####################################################################
# Replication code for
#   Hirose, Imai, and Lyall, "Can civilian attitudes predict insurgent violence?"
#
# Code to evaluate out-of-sample prediction performance 
#
# Author: Kentaro Hirose
####################################################################





########################################################################
library(endorse)
library(MASS)
library(mgcv)
library(sandwich)
library(robust)
library(Hmisc)
library(MCMCpack)
########################################################################
load("village.data.RData") ### in-sample village data
########################################################################
load("village.data.outsample.RData") ### out-sample village data
########################################################################
### in-sample violence (against ISAF)
load("pre.ied.attack.RData")
load("pre.ied.found.RData")
load("pre.nonied.RData")
load("post.ied.attack.RData")
load("post.ied.found.RData")
load("post.nonied.RData")
########################################################################
### in-sample violence (against civilians)
load("pre.ied.attack.civilian.RData")
load("pre.nonied.civilian.RData")
load("post.ied.attack.civilian.RData")
load("post.nonied.civilian.RData")
########################################################################
### out-sample violence (against ISAF)
load("other.pre.ied.attack.RData")
load("other.pre.ied.found.RData")
load("other.pre.nonied.RData")
load("other.post.ied.attack.RData")
load("other.post.ied.found.RData")
load("other.post.nonied.RData")
########################################################################
### out-sample violence (against civilians)
load("other.pre.ied.attack.civilian.RData")
load("other.pre.nonied.civilian.RData")
load("other.post.ied.attack.civilian.RData")
load("other.post.nonied.civilian.RData")
########################################################################








########################################################################
########################################################################
### Figure 5: out-of-sample prediction performance
########################################################################
########################################################################
verbose <- 1
K <- 60
T <- 10
MONTH <- 1:10
########################################################################
rate.taliban.ied.attack <- matrix(NA, T, K)
rate.taliban.ied.found <- matrix(NA, T, K)
rate.taliban.nonied <- matrix(NA, T, K)

rate.taliban.ied.attack.2 <- matrix(NA, T, K)
rate.taliban.ied.found.2 <- matrix(NA, T, K)
rate.taliban.nonied.2 <- matrix(NA, T, K)

ewA1 <- ewB1 <- ewC1 <- ewoA1 <- ewoB1 <- ewoC1 <- ewAA1 <- ewBB1 <- ewCC1 <- matrix(NA, T, K)
########################################################################
for (k in 1:K){
  for (t in 1:T){
    taliban <- village.data$support.taliban
    isaf <- village.data$support.isaf
    pre.A <- pre.ied.attack[[t]][, k]
    pre.B <- pre.ied.found[[t]][, k]
    pre.C <- pre.nonied[[t]][, k]
    post.A <- post.ied.attack[[t]][, k]
    post.B <- post.ied.found[[t]][, k]
    post.C <- post.nonied[[t]][, k]
    data <- data.frame(pre.A, post.A, pre.B, post.B, pre.C, post.C, taliban, isaf, village.data)
    IN <- na.omit(data)
    
    model <- lm(taliban ~
                village.data$log.population +
                village.data$log.altitude +
                village.data$control + 
                village.data$pakistan +
                village.data$talibsharia06)
    beta.taliban <- coef(model)

    model <- lm(isaf ~
                village.data$log.population +
                village.data$log.altitude +
                village.data$control +
                village.data$pakistan +
                village.data$talibsharia06)
    beta.isaf <- coef(model)


    pre.A <- other.pre.ied.attack[[t]][, k]
    pre.B <- other.pre.ied.found[[t]][, k]
    pre.C <- other.pre.nonied[[t]][, k]
    post.A <- other.post.ied.attack[[t]][, k]
    post.B <- other.post.ied.found[[t]][, k]
    post.C <- other.post.nonied[[t]][, k]
    
    X <- cbind(1,
               log(out.PPL$population),
               log(out.PPL$altitude),
               out.PPL$control,
               out.PPL$pakistan,
               out.PPL$talibsharia06)
    
    isaf <- as.vector(beta.isaf %*% t(X))
    taliban <- as.vector(beta.taliban %*% t(X))
    
    data <- data.frame(pre.A, post.A, pre.B, post.B, pre.C, post.C, taliban, isaf, out.PPL)
    OUT <- data
    OUT <- OUT[out.PPL$population != 0, ]
########################################################################
    w <- wo <- rep(NA, nrow(OUT))
########################################################################
    Y <- IN$post.A
    Y.lag <- IN$pre.A
    XX <- X<- IN$isaf - IN$taliban
    aid <- IN$aid
    base <- IN$base_3km
    out <- lm(Y ~ X + Y.lag + base + aid)
    w.beta <- coef(out)
    out <- lm(Y ~ Y.lag + base + aid)
    wo.beta <- coef(out)
    YY <- Y <- OUT$post.A
    YY.lag <- Y.lag <- OUT$pre.A
    X<- OUT$isaf - OUT$taliban   
    X <- X * (sd(XX)/sd(X))###
    aid <- OUT$aid
    base <- OUT$base_3km
    Z <- cbind(1, X, Y.lag, base, aid)
    w.mu <- as.vector(w.beta %*% t(Z))
    Z <- cbind(1, Y.lag, base, aid)
    wo.mu <- as.vector(wo.beta %*% t(Z))
    ewA1[t, k] <- w <- mean(abs(Y - w.mu), na.rm = TRUE)
    ewoA1[t, k] <- wo <- mean(abs(Y - wo.mu), na.rm = TRUE)  
    A1 <- rate.taliban.ied.attack[t, k] <- (wo-w)*100/w

    Y <- IN$post.B 
    Y.lag <- IN$pre.B 
    XX <- X<- IN$isaf - IN$taliban
    aid <- IN$aid
    base <- IN$base_3km
    out <- lm(Y ~ X + Y.lag + base + aid)
    w.beta <- coef(out)
    out <- lm(Y ~ Y.lag + base + aid)
    wo.beta <- coef(out)
    YY <- Y <- OUT$post.B 
    YY.lag <- Y.lag <- OUT$pre.B 
    X<- OUT$isaf - OUT$taliban   
    X <- X * (sd(XX)/sd(X))###
    aid <- OUT$aid
    base <- OUT$base_3km
    Z <- cbind(1, X, Y.lag, base, aid)
    w.mu <- as.vector(w.beta %*% t(Z))
    Z <- cbind(1, Y.lag, base, aid)
    wo.mu <- as.vector(wo.beta %*% t(Z))
    ewB1[t, k] <- w <- mean(abs(Y - w.mu), na.rm = TRUE)
    ewoB1[t, k] <- wo <- mean(abs(Y - wo.mu), na.rm = TRUE) 
    B1 <- rate.taliban.ied.found[t, k] <- (wo-w)*100/w

    Y <- IN$post.C 
    Y.lag <- IN$pre.C 
    XX <- X<- IN$isaf - IN$taliban
    aid <- IN$aid
    base <- IN$base_3km
    out <- lm(Y ~ X + Y.lag + base + aid)
    w.beta <- coef(out)
    out <- lm(Y ~ Y.lag + base + aid)
    wo.beta <- coef(out)
    YY <- Y <- OUT$post.C 
    YY.lag <- Y.lag <- OUT$pre.C 
    X<- OUT$isaf - OUT$taliban   
    X <- X * (sd(XX)/sd(X))###
    aid <- OUT$aid
    base <- OUT$base_3km
    Z <- cbind(1, X, Y.lag, base, aid)
    w.mu <- as.vector(w.beta %*% t(Z))
    Z <- cbind(1, Y.lag, base, aid)
    wo.mu <- as.vector(wo.beta %*% t(Z))
    ewC1[t, k] <- w <- mean(abs(Y - w.mu), na.rm = TRUE)
    ewoC1[t, k] <- wo <- mean(abs(Y - wo.mu), na.rm = TRUE) 
    C1 <- rate.taliban.nonied[t, k] <- (wo-w)*100/w

    
    Y <- IN$post.A
    Y.lag <- IN$pre.A
    aid <- IN$aid
    base <- IN$base_3km
    X<- cbind(IN$log.population, IN$log.altitude, IN$control, IN$pakistan, IN$talibsharia06) 
    out <- lm(Y ~ X + Y.lag + base + aid)
    w.beta <- coef(out)
    out <- lm(Y ~ Y.lag + base + aid)
    wo.beta <- coef(out)
    YY <- Y <- OUT$post.A
    YY.lag <- Y.lag <- OUT$pre.A
    aid <- OUT$aid
    base <- OUT$base_3km
    X<- cbind(log(OUT$population), log(OUT$altitude), OUT$control, OUT$pakistan, OUT$talibsharia06) 
    Z <- cbind(1, X, Y.lag, base, aid)
    w.mu <- as.vector(w.beta %*% t(Z))
    Z <- cbind(1, Y.lag, base, aid)
    wo.mu <- as.vector(wo.beta %*% t(Z))
    ewAA1[t, k] <- w <- mean(abs(Y - w.mu), na.rm = TRUE)
    wo <- mean(abs(Y - wo.mu), na.rm = TRUE) 
    A2 <- rate.taliban.ied.attack.2[t, k] <- (wo-w)*100/w

    Y <- IN$post.B 
    Y.lag <- IN$pre.B 
    aid <- IN$aid
    base <- IN$base_3km
    X<- cbind(IN$log.population, IN$log.altitude, IN$control, IN$pakistan, IN$talibsharia06) 
    out <- lm(Y ~ X + Y.lag + base + aid)
    w.beta <- coef(out)
    out <- lm(Y ~ Y.lag + base + aid)
    wo.beta <- coef(out)
    YY <- Y <- OUT$post.B 
    YY.lag <- Y.lag <- OUT$pre.B 
    aid <- OUT$aid
    base <- OUT$base_3km
    X<- cbind(log(OUT$population), log(OUT$altitude), OUT$control, OUT$pakistan, OUT$talibsharia06)
    Z <- cbind(1, X, Y.lag, base, aid)
    w.mu <- as.vector(w.beta %*% t(Z))
    Z <- cbind(1, Y.lag, base, aid)
    wo.mu <- as.vector(wo.beta %*% t(Z))
    ewBB1[t, k] <- w <- mean(abs(Y - w.mu), na.rm = TRUE)
    wo <- mean(abs(Y - wo.mu), na.rm = TRUE) 
    B2 <- rate.taliban.ied.found.2[t, k] <- (wo-w)*100/w

    Y <- IN$post.C 
    Y.lag <- IN$pre.C 
    aid <- IN$aid
    base <- IN$base_3km
    X<- cbind(IN$log.population, IN$log.altitude, IN$control, IN$pakistan, IN$talibsharia06) 
    out <- lm(Y ~ X + Y.lag + base + aid)
    w.beta <- coef(out)
    out <- lm(Y ~ Y.lag + base + aid)
    wo.beta <- coef(out)
    YY <- Y <- OUT$post.C 
    YY.lag <- Y.lag <- OUT$pre.C 
    aid <- OUT$aid
    base <- OUT$base_3km
    X<- cbind(log(OUT$population), log(OUT$altitude), OUT$control, OUT$pakistan, OUT$talibsharia06)
    Z <- cbind(1, X, Y.lag, base, aid)
    w.mu <- as.vector(w.beta %*% t(Z))
    Z <- cbind(1, Y.lag, base, aid)
    wo.mu <- as.vector(wo.beta %*% t(Z))
    ewCC1[t, k] <- w <- mean(abs(Y - w.mu), na.rm = TRUE)
    wo <- mean(abs(Y - wo.mu), na.rm = TRUE) 
    C2 <- rate.taliban.nonied.2[t, k] <- (wo-w)*100/w
    
    
    if (verbose!=0&t%%verbose == 0){
      cat("t", t, "k", k, "A1", A1, "B1", B1, "C1", C1, "A2", A2, "B2", B2, "C2", C2, "\n") 
    }
  }
}
##############################################################
png(file = "pred_support_vs_covariate.png", width = 1800, height = 1500)


ccc <- 3
mmm <- 3
ddd <- 3
fff <- 3
XLAB <- list("Distance window (km)", cex = ccc)
YLAB <- list("Time window (months)", cex = ccc)
COL.MIN <- -43
COL.MAX <- 43
BY <- 2
BY2 <- 2

Q <- rate.taliban.ied.attack
MAIN <- list("IED attacks", cex = mmm)
elevation.df <- data.frame(x = rep(1:K, each = T), y = rep(seq(1, MONTH[T], 1), K), z = as.vector(Q))
elevation.loess <- loess(z ~ x * y, data = elevation.df)
elevation.fit <- expand.grid(list(x = seq(1, K, 0.05), y = seq(1, MONTH[T], 0.05)))
z <- predict(elevation.loess, newdata = elevation.fit)
elevation.fit$z <- as.numeric(z)
P1 <- levelplot(z ~ x*y, data = elevation.fit,
                xlab = XLAB, ylab = YLAB,
                main = MAIN,
                col.regions = colorRampPalette(c('dark red','white','dark blue')),
                contour = TRUE,
                at = seq(COL.MIN, COL.MAX, BY), 
                labels = list(TRUE, cex = fff),
                colorkey = FALSE,
                aspect = "xy", scales=list(x=list(cex = ddd),y=list(cex = ddd))
                )


Q <- rate.taliban.ied.found
MAIN <- list("IED found", cex = mmm)
elevation.df <- data.frame(x = rep(1:K, each = T), y = rep(seq(1, MONTH[T], 1), K), z = as.vector(Q))
elevation.loess <- loess(z ~ x * y, data = elevation.df)
elevation.fit <- expand.grid(list(x = seq(1, K, 0.05), y = seq(1, MONTH[T], 0.05)))
z <- predict(elevation.loess, newdata = elevation.fit)
elevation.fit$z <- as.numeric(z)
P2 <- levelplot(z ~ x*y, data = elevation.fit,
                xlab = XLAB, ylab = YLAB,
                main = MAIN,
                col.regions = colorRampPalette(c('dark red','white','dark blue')),
                contour = TRUE,
                at = seq(COL.MIN, COL.MAX, BY), 
                labels = list(TRUE, cex = fff),
                colorkey = FALSE,
                aspect = "xy", scales=list(x=list(cex = ddd),y=list(cex = ddd))
                )

Q <- rate.taliban.nonied
MAIN <- list("Non-IED attacks", cex = mmm)
elevation.df <- data.frame(x = rep(1:K, each = T), y = rep(seq(1, MONTH[T], 1), K), z = as.vector(Q))
elevation.loess <- loess(z ~ x * y, data = elevation.df)
elevation.fit <- expand.grid(list(x = seq(1, K, 0.05), y = seq(1, MONTH[T], 0.05)))
z <- predict(elevation.loess, newdata = elevation.fit)
elevation.fit$z <- as.numeric(z)
P3 <- levelplot(z ~ x*y, data = elevation.fit,
                xlab = XLAB, ylab = YLAB,
                main = MAIN,
                col.regions = colorRampPalette(c('dark red','white','dark blue')),
                contour = TRUE,
                at = seq(COL.MIN, COL.MAX, BY), 
                labels = list(TRUE, cex = fff),
                colorkey = list(TRUE, labels = list(cex = fff)),
                aspect = "xy", scales=list(x=list(cex = ddd),y=list(cex = ddd))
                )


BY <- 5
BY2 <- 5
Q <- rate.taliban.ied.attack.2
MAIN <- list("IED attacks", cex = mmm)
elevation.df <- data.frame(x = rep(1:K, each = T), y = rep(seq(1, MONTH[T], 1), K), z = as.vector(Q))
elevation.loess <- loess(z ~ x * y, data = elevation.df)
elevation.fit <- expand.grid(list(x = seq(1, K, 0.05), y = seq(1, MONTH[T], 0.05)))
z <- predict(elevation.loess, newdata = elevation.fit)
elevation.fit$z <- as.numeric(z)
P4 <- levelplot(z ~ x*y, data = elevation.fit,
                xlab = XLAB, ylab = YLAB,
                main = MAIN,
                col.regions = colorRampPalette(c('dark red','white','dark blue')),
                contour = TRUE,
                at = seq(COL.MIN, COL.MAX, BY), 
                labels = list(TRUE, cex = fff),
                colorkey = FALSE,
                aspect = "xy", scales=list(x=list(cex = ddd),y=list(cex = ddd))
                )


Q <- rate.taliban.ied.found.2
MAIN <- list("IED found", cex = mmm)
elevation.df <- data.frame(x = rep(1:K, each = T), y = rep(seq(1, MONTH[T], 1), K), z = as.vector(Q))
elevation.loess <- loess(z ~ x * y, data = elevation.df)
elevation.fit <- expand.grid(list(x = seq(1, K, 0.05), y = seq(1, MONTH[T], 0.05)))
z <- predict(elevation.loess, newdata = elevation.fit)
elevation.fit$z <- as.numeric(z)
P5 <- levelplot(z ~ x*y, data = elevation.fit,
                xlab = XLAB, ylab = YLAB,
                main = MAIN,
                col.regions = colorRampPalette(c('dark red','white','dark blue')),
                contour = TRUE,
                at = seq(COL.MIN, COL.MAX, BY), 
                labels = list(TRUE, cex = fff),
                colorkey = FALSE,
                aspect = "xy", scales=list(x=list(cex = ddd),y=list(cex = ddd))
                )

Q <- rate.taliban.nonied.2
MAIN <- list("Non-IED attacks", cex = mmm)
elevation.df <- data.frame(x = rep(1:K, each = T), y = rep(seq(1, MONTH[T], 1), K), z = as.vector(Q))
elevation.loess <- loess(z ~ x * y, data = elevation.df)
elevation.fit <- expand.grid(list(x = seq(1, K, 0.05), y = seq(1, MONTH[T], 0.05)))
z <- predict(elevation.loess, newdata = elevation.fit)
elevation.fit$z <- as.numeric(z)
P6 <- levelplot(z ~ x*y, data = elevation.fit,
                xlab = XLAB, ylab = YLAB,
                main = MAIN,
                col.regions = colorRampPalette(c('dark red','white','dark blue')),
                contour = TRUE,
                at = seq(COL.MIN, COL.MAX, BY), 
                labels = list(TRUE, cex = fff),
                colorkey = list(TRUE, labels = list(cex = fff)),
                aspect = "xy", scales=list(x=list(cex = ddd),y=list(cex = ddd))
                )

ad1 <- 0.02
print(P1, position = c(0, 0.5 + ad1, 0.315, 0.95), more = TRUE)
print(P2, position = c(0.315, 0.5 + ad1, 0.63, 0.95), more = TRUE)
print(P3, position = c(0.63, 0.5 + ad1, 1, 0.95), more = TRUE)
print(P4, position = c(0, 0 + ad1, 0.315, 0.45), more = TRUE)
print(P5, position = c(0.315, 0 + ad1, 0.63, 0.45), more = TRUE)
print(P6, position = c(0.63, 0 + ad1, 1, 0.45))

panel.text(x = 850, y = 40, labels = "Prediction improvement due to support measure", cex = 4, font = 1)
panel.text(x = 850, y = 800, labels = "Prediction improvement due to covariates", cex = 4, font = 1)


dev.off()
##############################################################










########################################################################
########################################################################
###  Figure 23: out-of-sample prediction performance (RMSE)
########################################################################
########################################################################
verbose <- 1
K <- 60
T <- 10
MONTH <- 1:10
########################################################################
rate.taliban.ied.attack <- matrix(NA, T, K)
rate.taliban.ied.found <- matrix(NA, T, K)
rate.taliban.nonied <- matrix(NA, T, K)

rate.taliban.ied.attack.2 <- matrix(NA, T, K)
rate.taliban.ied.found.2 <- matrix(NA, T, K)
rate.taliban.nonied.2 <- matrix(NA, T, K)

ewA <- ewB <- ewC <- ewoA <- ewoB <- ewoC <- ewAA <- ewBB <- ewCC <- matrix(NA, T, K)
########################################################################
for (k in 1:K){
  for (t in 1:T){
    taliban <- village.data$support.taliban
    isaf <- village.data$support.isaf
    pre.A <- pre.ied.attack[[t]][, k]
    pre.B <- pre.ied.found[[t]][, k]
    pre.C <- pre.nonied[[t]][, k]
    post.A <- post.ied.attack[[t]][, k]
    post.B <- post.ied.found[[t]][, k]
    post.C <- post.nonied[[t]][, k]
    data <- data.frame(pre.A, post.A, pre.B, post.B, pre.C, post.C, taliban, isaf, village.data)
    IN <- na.omit(data)
    
    model <- lm(taliban ~
                village.data$log.population +
                village.data$log.altitude +
                village.data$control +
                village.data$pakistan +
                village.data$talibsharia06)
    beta.taliban <- coef(model)

    model <- lm(isaf ~
                village.data$log.population +
                village.data$log.altitude +
                village.data$control +
                village.data$pakistan +
                village.data$talibsharia06)
    beta.isaf <- coef(model)


    pre.A <- other.pre.ied.attack[[t]][, k]
    pre.B <- other.pre.ied.found[[t]][, k]
    pre.C <- other.pre.nonied[[t]][, k]
    post.A <- other.post.ied.attack[[t]][, k]
    post.B <- other.post.ied.found[[t]][, k]
    post.C <- other.post.nonied[[t]][, k]
    
    X <- cbind(1,
               log(out.PPL$population),
               log(out.PPL$altitude),
               out.PPL$control,
               out.PPL$pakistan,
               out.PPL$talibsharia06)
    
    isaf <- as.vector(beta.isaf %*% t(X))
    taliban <- as.vector(beta.taliban %*% t(X))
    
    data <- data.frame(pre.A, post.A, pre.B, post.B, pre.C, post.C, taliban, isaf, out.PPL)
    OUT <- data
    OUT <- OUT[out.PPL$population != 0, ]
########################################################################
### With Support
########################################################################
    w <- wo <- rep(NA, nrow(OUT))
########################################################################
    Y <- IN$post.A
    Y.lag <- IN$pre.A
    XX <- X<- IN$isaf - IN$taliban
    aid <- IN$aid
    base <- IN$base_3km
    out <- lm(Y ~ X + Y.lag + base + aid)
    w.beta <- coef(out)
    out <- lm(Y ~ Y.lag + base + aid)
    wo.beta <- coef(out)
    YY <- Y <- OUT$post.A
    YY.lag <- Y.lag <- OUT$pre.A
    X<- OUT$isaf - OUT$taliban   
    X <- X * (sd(XX)/sd(X))###
    aid <- OUT$aid
    base <- OUT$base_3km
    Z <- cbind(1, X, Y.lag, base, aid)
    w.mu <- as.vector(w.beta %*% t(Z))
    Z <- cbind(1, Y.lag, base, aid)
    wo.mu <- as.vector(wo.beta %*% t(Z))
    w <- sqrt(mean((Y - w.mu)^2, na.rm = TRUE))
    wo <- sqrt(mean((Y - wo.mu)^2, na.rm = TRUE)) 
    A1 <- rate.taliban.ied.attack[t, k] <- (wo-w)*100/w

    Y <- IN$post.B 
    Y.lag <- IN$pre.B 
    XX <- X<- IN$isaf - IN$taliban
    aid <- IN$aid
    base <- IN$base_3km
    out <- lm(Y ~ X + Y.lag + base + aid)
    w.beta <- coef(out)
    out <- lm(Y ~ Y.lag + base + aid)
    wo.beta <- coef(out)
    YY <- Y <- OUT$post.B 
    YY.lag <- Y.lag <- OUT$pre.B 
    X<- OUT$isaf - OUT$taliban   
    X <- X * (sd(XX)/sd(X))###
    aid <- OUT$aid
    base <- OUT$base_3km
    Z <- cbind(1, X, Y.lag, base, aid)
    w.mu <- as.vector(w.beta %*% t(Z))
    Z <- cbind(1, Y.lag, base, aid)
    wo.mu <- as.vector(wo.beta %*% t(Z))
    w <- sqrt(mean((Y - w.mu)^2, na.rm = TRUE))
    wo <- sqrt(mean((Y - wo.mu)^2, na.rm = TRUE)) 
    B1 <- rate.taliban.ied.found[t, k] <- (wo-w)*100/w

    Y <- IN$post.C 
    Y.lag <- IN$pre.C 
    XX <- X<- IN$isaf - IN$taliban
    aid <- IN$aid
    base <- IN$base_3km
    out <- lm(Y ~ X + Y.lag + base + aid)
    w.beta <- coef(out)
    out <- lm(Y ~ Y.lag + base + aid)
    wo.beta <- coef(out)
    YY <- Y <- OUT$post.C 
    YY.lag <- Y.lag <- OUT$pre.C 
    X<- OUT$isaf - OUT$taliban   
    X <- X * (sd(XX)/sd(X))###
    aid <- OUT$aid
    base <- OUT$base_3km
    Z <- cbind(1, X, Y.lag, base, aid)
    w.mu <- as.vector(w.beta %*% t(Z))
    Z <- cbind(1, Y.lag, base, aid)
    wo.mu <- as.vector(wo.beta %*% t(Z))
    w <- sqrt(mean((Y - w.mu)^2, na.rm = TRUE))
    wo <- sqrt(mean((Y - wo.mu)^2, na.rm = TRUE)) 
    C1 <- rate.taliban.nonied[t, k] <- (wo-w)*100/w

    
    Y <- IN$post.A
    Y.lag <- IN$pre.A
    aid <- IN$aid
    base <- IN$base_3km
    X<- cbind(IN$log.population, IN$log.altitude, IN$control, IN$pakistan, IN$talibsharia06) 
    out <- lm(Y ~ X + Y.lag + base + aid)
    w.beta <- coef(out)
    out <- lm(Y ~ Y.lag + base + aid)
    wo.beta <- coef(out)
    YY <- Y <- OUT$post.A
    YY.lag <- Y.lag <- OUT$pre.A
    aid <- OUT$aid
    base <- OUT$base_3km
    X<- cbind(log(OUT$population), log(OUT$altitude), OUT$control, OUT$pakistan, OUT$talibsharia06)  
    Z <- cbind(1, X, Y.lag, base, aid)
    w.mu <- as.vector(w.beta %*% t(Z))
    Z <- cbind(1, Y.lag, base, aid)
    wo.mu <- as.vector(wo.beta %*% t(Z))
    w <- sqrt(mean((Y - w.mu)^2, na.rm = TRUE))
    wo <- sqrt(mean((Y - wo.mu)^2, na.rm = TRUE)) 
    A2 <- rate.taliban.ied.attack.2[t, k] <- (wo-w)*100/w

    Y <- IN$post.B 
    Y.lag <- IN$pre.B 
    aid <- IN$aid
    base <- IN$base_3km
    X<- cbind(IN$log.population, IN$log.altitude, IN$control, IN$pakistan, IN$talibsharia06) 
    out <- lm(Y ~ X + Y.lag + base + aid)
    w.beta <- coef(out)
    out <- lm(Y ~ Y.lag + base + aid)
    wo.beta <- coef(out)
    YY <- Y <- OUT$post.B 
    YY.lag <- Y.lag <- OUT$pre.B 
    aid <- OUT$aid
    base <- OUT$base_3km
    X<- cbind(log(OUT$population), log(OUT$altitude), OUT$control, OUT$pakistan, OUT$talibsharia06)
    Z <- cbind(1, X, Y.lag, base, aid)
    w.mu <- as.vector(w.beta %*% t(Z))
    Z <- cbind(1, Y.lag, base, aid)
    wo.mu <- as.vector(wo.beta %*% t(Z))
    w <- sqrt(mean((Y - w.mu)^2, na.rm = TRUE))
    wo <- sqrt(mean((Y - wo.mu)^2, na.rm = TRUE)) 
    B2 <- rate.taliban.ied.found.2[t, k] <- (wo-w)*100/w

    Y <- IN$post.C 
    Y.lag <- IN$pre.C 
    aid <- IN$aid
    base <- IN$base_3km
    X<- cbind(IN$log.population, IN$log.altitude, IN$control, IN$pakistan, IN$talibsharia06) 
    out <- lm(Y ~ X + Y.lag + base + aid)
    w.beta <- coef(out)
    out <- lm(Y ~ Y.lag + base + aid)
    wo.beta <- coef(out)
    YY <- Y <- OUT$post.C 
    YY.lag <- Y.lag <- OUT$pre.C 
    aid <- OUT$aid
    base <- OUT$base_3km
    X<- cbind(log(OUT$population), log(OUT$altitude), OUT$control, OUT$pakistan, OUT$talibsharia06)
    Z <- cbind(1, X, Y.lag, base, aid)
    w.mu <- as.vector(w.beta %*% t(Z))
    Z <- cbind(1, Y.lag, base, aid)
    wo.mu <- as.vector(wo.beta %*% t(Z))
    w <- sqrt(mean((Y - w.mu)^2, na.rm = TRUE))
    wo <- sqrt(mean((Y - wo.mu)^2, na.rm = TRUE)) 
    C2 <- rate.taliban.nonied.2[t, k] <- (wo-w)*100/w
    
    
    if (verbose!=0&t%%verbose == 0){
      cat("t", t, "k", k, "A1", A1, "B1", B1, "C1", C1, "A2", A2, "B2", B2, "C2", C2, "\n") 
    }
  }
}
##############################################################
png(file = "pred_rmse_support_vs_covariate.png", width = 1800, height =1500)


ccc <- 3
mmm <- 3
ddd <- 3
fff <- 3
XLAB <- list("Distance window (km)", cex = ccc)
YLAB <- list("Time window (months)", cex = ccc)
COL.MIN <- -30
COL.MAX <- 30
BY <- 2
BY2 <- 2

Q <- rate.taliban.ied.attack
MAIN <- list("IED attacks", cex = mmm)
elevation.df <- data.frame(x = rep(1:K, each = T), y = rep(seq(1, MONTH[T], 1), K), z = as.vector(Q))
elevation.loess <- loess(z ~ x * y, data = elevation.df)
elevation.fit <- expand.grid(list(x = seq(1, K, 0.05), y = seq(1, MONTH[T], 0.05)))
z <- predict(elevation.loess, newdata = elevation.fit)
elevation.fit$z <- as.numeric(z)
P1 <- levelplot(z ~ x*y, data = elevation.fit,
                xlab = XLAB, ylab = YLAB,
                main = MAIN,
                col.regions = colorRampPalette(c('dark red','white','dark blue')),
                contour = TRUE,
                at = seq(COL.MIN, COL.MAX, BY), 
                labels = list(TRUE, cex = fff),
                colorkey = FALSE,
                aspect = "xy", scales=list(x=list(cex = ddd),y=list(cex = ddd))
                )

Q <- rate.taliban.ied.found
MAIN <- list("IED found", cex = mmm)
elevation.df <- data.frame(x = rep(1:K, each = T), y = rep(seq(1, MONTH[T], 1), K), z = as.vector(Q))
elevation.loess <- loess(z ~ x * y, data = elevation.df)
elevation.fit <- expand.grid(list(x = seq(1, K, 0.05), y = seq(1, MONTH[T], 0.05)))
z <- predict(elevation.loess, newdata = elevation.fit)
elevation.fit$z <- as.numeric(z)
P2 <- levelplot(z ~ x*y, data = elevation.fit,
                xlab = XLAB, ylab = YLAB,
                main = MAIN,
                col.regions = colorRampPalette(c('dark red','white','dark blue')),
                contour = TRUE,
                at = seq(COL.MIN, COL.MAX, BY), 
                labels = list(TRUE, cex = fff),
                colorkey = FALSE,
                aspect = "xy", scales=list(x=list(cex = ddd),y=list(cex = ddd))
                )

Q <- rate.taliban.nonied
MAIN <- list("Non-IED attacks", cex = mmm)
elevation.df <- data.frame(x = rep(1:K, each = T), y = rep(seq(1, MONTH[T], 1), K), z = as.vector(Q))
elevation.loess <- loess(z ~ x * y, data = elevation.df)
elevation.fit <- expand.grid(list(x = seq(1, K, 0.05), y = seq(1, MONTH[T], 0.05)))
z <- predict(elevation.loess, newdata = elevation.fit)
elevation.fit$z <- as.numeric(z)
P3 <- levelplot(z ~ x*y, data = elevation.fit,
                xlab = XLAB, ylab = YLAB,
                main = MAIN,
                col.regions = colorRampPalette(c('dark red','white','dark blue')),
                contour = TRUE,
                at = seq(COL.MIN, COL.MAX, BY), 
                labels = list(TRUE, cex = fff),
                colorkey = list(TRUE, labels = list(cex = fff)),
                aspect = "xy", scales=list(x=list(cex = ddd),y=list(cex = ddd))
                )

BY <- 5
BY2 <- 5
Q <- rate.taliban.ied.attack.2
MAIN <- list("IED attacks", cex = mmm)
elevation.df <- data.frame(x = rep(1:K, each = T), y = rep(seq(1, MONTH[T], 1), K), z = as.vector(Q))
elevation.loess <- loess(z ~ x * y, data = elevation.df)
elevation.fit <- expand.grid(list(x = seq(1, K, 0.05), y = seq(1, MONTH[T], 0.05)))
z <- predict(elevation.loess, newdata = elevation.fit)
elevation.fit$z <- as.numeric(z)
P4 <- levelplot(z ~ x*y, data = elevation.fit,
                xlab = XLAB, ylab = YLAB,
                main = MAIN,
                col.regions = colorRampPalette(c('dark red','white','dark blue')),
                contour = TRUE,
                at = seq(COL.MIN, COL.MAX, BY), 
                labels = list(TRUE, cex = fff),
                colorkey = FALSE,
                aspect = "xy", scales=list(x=list(cex = ddd),y=list(cex = ddd))
                )


Q <- rate.taliban.ied.found.2
MAIN <- list("IED found", cex = mmm)
elevation.df <- data.frame(x = rep(1:K, each = T), y = rep(seq(1, MONTH[T], 1), K), z = as.vector(Q))
elevation.loess <- loess(z ~ x * y, data = elevation.df)
elevation.fit <- expand.grid(list(x = seq(1, K, 0.05), y = seq(1, MONTH[T], 0.05)))
z <- predict(elevation.loess, newdata = elevation.fit)
elevation.fit$z <- as.numeric(z)
P5 <- levelplot(z ~ x*y, data = elevation.fit,
                xlab = XLAB, ylab = YLAB,
                main = MAIN,
                col.regions = colorRampPalette(c('dark red','white','dark blue')),
                contour = TRUE,
                at = seq(COL.MIN, COL.MAX, BY), 
                labels = list(TRUE, cex = fff),
                colorkey = FALSE,
                aspect = "xy", scales=list(x=list(cex = ddd),y=list(cex = ddd))
                )

Q <- rate.taliban.nonied.2
MAIN <- list("Non-IED attacks", cex = mmm)
elevation.df <- data.frame(x = rep(1:K, each = T), y = rep(seq(1, MONTH[T], 1), K), z = as.vector(Q))
elevation.loess <- loess(z ~ x * y, data = elevation.df)
elevation.fit <- expand.grid(list(x = seq(1, K, 0.05), y = seq(1, MONTH[T], 0.05)))
z <- predict(elevation.loess, newdata = elevation.fit)
elevation.fit$z <- as.numeric(z)
P6 <- levelplot(z ~ x*y, data = elevation.fit,
                xlab = XLAB, ylab = YLAB,
                main = MAIN,
                col.regions = colorRampPalette(c('dark red','white','dark blue')),
                contour = TRUE,
                at = seq(COL.MIN, COL.MAX, BY), 
                labels = list(TRUE, cex = fff),
                colorkey = list(TRUE, labels = list(cex = fff)),
                aspect = "xy", scales=list(x=list(cex = ddd),y=list(cex = ddd))
                )

ad1 <- 0.02
print(P1, position = c(0, 0.5 + ad1, 0.315, 0.95), more = TRUE)
print(P2, position = c(0.315, 0.5 + ad1, 0.63, 0.95), more = TRUE)
print(P3, position = c(0.63, 0.5 + ad1, 1, 0.95), more = TRUE)
print(P4, position = c(0, 0 + ad1, 0.315, 0.45), more = TRUE)
print(P5, position = c(0.315, 0 + ad1, 0.63, 0.45), more = TRUE)
print(P6, position = c(0.63, 0 + ad1, 1, 0.45))

panel.text(x = 850, y = 40, labels = "Prediction improvement due to support measure", cex = 4, font = 1)
panel.text(x = 850, y = 800, labels = "Prediction improvement due to covariates", cex = 4, font = 1)


dev.off()
##############################################################




########################################################################
########################################################################
###  Figure 24: Interaction term I
########################################################################
########################################################################
verbose <- 1
K <- 60
T <- 10
MONTH <- 1:10
########################################################################
rate.taliban.ied.attack.1 <- matrix(NA, T, K)
rate.taliban.ied.found.1 <- matrix(NA, T, K)
rate.taliban.nonied.1 <- matrix(NA, T, K)

rate.taliban.ied.attack.2 <- matrix(NA, T, K)
rate.taliban.ied.found.2 <- matrix(NA, T, K)
rate.taliban.nonied.2 <- matrix(NA, T, K)

rate.taliban.ied.attack.3 <- matrix(NA, T, K)
rate.taliban.ied.found.3 <- matrix(NA, T, K)
rate.taliban.nonied.3 <- matrix(NA, T, K)
########################################################################
for (k in 1:K){
  for (t in 1:T){
    taliban <- village.data$support.taliban
    isaf <- village.data$support.isaf
    pre.A <- pre.ied.attack[[t]][, k]
    pre.B <- pre.ied.found[[t]][, k]
    pre.C <- pre.nonied[[t]][, k]
    post.A <- post.ied.attack[[t]][, k]
    post.B <- post.ied.found[[t]][, k]
    post.C <- post.nonied[[t]][, k]
    data <- data.frame(pre.A, post.A, pre.B, post.B, pre.C, post.C, taliban, isaf, village.data)
    IN <- na.omit(data)
    
    model <- lm(taliban ~
                village.data$log.population +
                village.data$log.altitude +
                village.data$control +
                village.data$pakistan +
                village.data$talibsharia06)
    beta.taliban <- coef(model)

    model <- lm(isaf ~
                village.data$log.population +
                village.data$log.altitude +
                village.data$control +
                village.data$pakistan +
                village.data$talibsharia06)
    beta.isaf <- coef(model)


    pre.A <- other.pre.ied.attack[[t]][, k]
    pre.B <- other.pre.ied.found[[t]][, k]
    pre.C <- other.pre.nonied[[t]][, k]
    post.A <- other.post.ied.attack[[t]][, k]
    post.B <- other.post.ied.found[[t]][, k]
    post.C <- other.post.nonied[[t]][, k]
    
    X <- cbind(1,
           log(out.PPL$population),
           log(out.PPL$altitude),
           out.PPL$control,
           out.PPL$pakistan,
           out.PPL$talibsharia06)
    
    isaf <- as.vector(beta.isaf %*% t(X))
    taliban <- as.vector(beta.taliban %*% t(X))
    
    data <- data.frame(pre.A, post.A, pre.B, post.B, pre.C, post.C, taliban, isaf, out.PPL)
    OUT <- data
    OUT <- OUT[out.PPL$population != 0, ]
########################################################################
    w <- wo <- rep(NA, nrow(OUT))
########################################################################
    Y <- IN$post.A
    Y.lag <- IN$pre.A
    XX <- X<- IN$isaf - IN$taliban
    aid <- IN$aid
    base <- IN$base_3km
    
    #out <- lm(Y ~ X + Y.lag + base + aid + X:base + X:aid)
    #out <- lm(Y ~ X + Y.lag + base + aid + X:Y.lag + X:aid)
    #out <- lm(Y ~ X + Y.lag + base + aid + X:Y.lag + X:base)
    out <- lm(Y ~ X + Y.lag + base + aid + X:Y.lag)
    #out <- lm(Y ~ X + Y.lag + base + aid + X:base)
    #out <- lm(Y ~ X + Y.lag + base + aid + X:aid)
    #out <- lm(Y ~ X + Y.lag + base + aid)
    
    w.beta <- coef(out)
    out <- lm(Y ~ Y.lag + base + aid)
    wo.beta <- coef(out)
    Y <- OUT$post.A
    Y.lag <- OUT$pre.A
    X<- OUT$isaf - OUT$taliban
    X <- X * (sd(XX)/sd(X))### 
    aid <- OUT$aid
    base <- OUT$base_3km

    #Z <- cbind(1, X, Y.lag, base, aid, X*base, X*aid)
    #Z <- cbind(1, X, Y.lag, base, aid, X*Y.lag, X*aid)
    #Z <- cbind(1, X, Y.lag, base, aid, X*Y.lag, X*base)
    Z <- cbind(1, X, Y.lag, base, aid, X*Y.lag)
    #Z <- cbind(1, X, Y.lag, base, aid, X*base)
    #Z <- cbind(1, X, Y.lag, base, aid, X*aid)
    #Z <- cbind(1, X, Y.lag, base, aid)    
    
    w.mu <- as.vector(w.beta %*% t(Z))
    Z <- cbind(1, Y.lag, base, aid)
    wo.mu <- as.vector(wo.beta %*% t(Z))
    w <- mean(abs(Y - w.mu))
    wo <- mean(abs(Y - wo.mu))  
    A1 <- rate.taliban.ied.attack.1[t, k] <- (wo-w)*100/w
    
    ####################################################
    Y <- IN$post.B 
    Y.lag <- IN$pre.B 
    XX <- X <- IN$isaf - IN$taliban
    aid <- IN$aid
    base <- IN$base_3km

    #out <- lm(Y ~ X + Y.lag + base + aid + X:base + X:aid)
    #out <- lm(Y ~ X + Y.lag + base + aid + X:Y.lag + X:aid)
    #out <- lm(Y ~ X + Y.lag + base + aid + X:Y.lag + X:base)
    out <- lm(Y ~ X + Y.lag + base + aid + X:Y.lag)
    #out <- lm(Y ~ X + Y.lag + base + aid + X:base)
    #out <- lm(Y ~ X + Y.lag + base + aid + X:aid)
    #out <- lm(Y ~ X + Y.lag + base + aid + X:aid)
    
    w.beta <- coef(out)
    out <- lm(Y ~ Y.lag + base + aid)
    wo.beta <- coef(out)
    Y <- OUT$post.B 
    Y.lag <- OUT$pre.B 
    X<- OUT$isaf - OUT$taliban
    X <- X * (sd(XX)/sd(X))###
    aid <- OUT$aid
    base <- OUT$base_3km
    
    #Z <- cbind(1, X, Y.lag, base, aid, X*base, X*aid)
    #Z <- cbind(1, X, Y.lag, base, aid, X*Y.lag, X*aid)
    #Z <- cbind(1, X, Y.lag, base, aid, X*Y.lag, X*base)
    Z <- cbind(1, X, Y.lag, base, aid, X*Y.lag)
    #Z <- cbind(1, X, Y.lag, base, aid, X*base)
    #Z <- cbind(1, X, Y.lag, base, aid, X*aid)
    #Z <- cbind(1, X, Y.lag, base, aid, X*aid)
    
    w.mu <- w.beta %*% t(Z)
    Z <- cbind(1, Y.lag, base, aid)
    wo.mu <- wo.beta %*% t(Z)
    w <- mean(abs(Y - w.mu))
    wo <- mean(abs(Y - wo.mu))  
    B1 <- rate.taliban.ied.found.1[t, k] <- (wo-w)*100/w

    ####################################################
    Y <- IN$post.C 
    Y.lag <- IN$pre.C 
    XX <- X <- IN$isaf - IN$taliban
    aid <- IN$aid
    base <- IN$base_3km

    #out <- lm(Y ~ X + Y.lag + base + aid + X:base + X:aid)
    #out <- lm(Y ~ X + Y.lag + base + aid + X:Y.lag + X:aid)
    #out <- lm(Y ~ X + Y.lag + base + aid + X:Y.lag + X:base)
    out <- lm(Y ~ X + Y.lag + base + aid + X:Y.lag)
    #out <- lm(Y ~ X + Y.lag + base + aid + X:base)
    #out <- lm(Y ~ X + Y.lag + base + aid + X:aid)
    #out <- lm(Y ~ X + Y.lag + base + aid + X:Y.lag)
    
    w.beta <- coef(out)
    out <- lm(Y ~ Y.lag + base + aid)
    wo.beta <- coef(out)
    Y <- OUT$post.C 
    Y.lag <- OUT$pre.C 
    X<- OUT$isaf - OUT$taliban
    X <- X * (sd(XX)/sd(X))###
    aid <- OUT$aid
    base <- OUT$base_3km
    
    #Z <- cbind(1, X, Y.lag, base, aid, X*base, X*aid)
    #Z <- cbind(1, X, Y.lag, base, aid, X*Y.lag, X*aid)
    #Z <- cbind(1, X, Y.lag, base, aid, X*Y.lag, X*base)
    Z <- cbind(1, X, Y.lag, base, aid, X*Y.lag)
    #Z <- cbind(1, X, Y.lag, base, aid, X*base)
    #Z <- cbind(1, X, Y.lag, base, aid, X*aid)
    #Z <- cbind(1, X, Y.lag, base, aid, X*Y.lag)
    
    w.mu <- w.beta %*% t(Z)
    Z <- cbind(1, Y.lag, base, aid)
    wo.mu <- wo.beta %*% t(Z)
    w <- mean(abs(Y - w.mu))
    wo <- mean(abs(Y - wo.mu))  
    C1 <- rate.taliban.nonied.1[t, k] <- (wo-w)*100/w


########################################################################
    Y <- IN$post.A
    Y.lag <- IN$pre.A
    XX <- X<- IN$isaf - IN$taliban
    aid <- IN$aid
    base <- IN$base_3km
    
    #out <- lm(Y ~ X + Y.lag + base + aid + X:base + X:aid)
    #out <- lm(Y ~ X + Y.lag + base + aid + X:Y.lag + X:aid)
    #out <- lm(Y ~ X + Y.lag + base + aid + X:Y.lag + X:base)
    #out <- lm(Y ~ X + Y.lag + base + aid + X:Y.lag)
    out <- lm(Y ~ X + Y.lag + base + aid + X:base)
    #out <- lm(Y ~ X + Y.lag + base + aid + X:aid)
    #out <- lm(Y ~ X + Y.lag + base + aid)
    
    w.beta <- coef(out)
    out <- lm(Y ~ Y.lag + base + aid)
    wo.beta <- coef(out)
    Y <- OUT$post.A
    Y.lag <- OUT$pre.A
    X<- OUT$isaf - OUT$taliban
    X <- X * (sd(XX)/sd(X))### 
    aid <- OUT$aid
    base <- OUT$base_3km

    #Z <- cbind(1, X, Y.lag, base, aid, X*base, X*aid)
    #Z <- cbind(1, X, Y.lag, base, aid, X*Y.lag, X*aid)
    #Z <- cbind(1, X, Y.lag, base, aid, X*Y.lag, X*base)
    #Z <- cbind(1, X, Y.lag, base, aid, X*Y.lag)
    Z <- cbind(1, X, Y.lag, base, aid, X*base)
    #Z <- cbind(1, X, Y.lag, base, aid, X*aid)
    #Z <- cbind(1, X, Y.lag, base, aid)    
    
    w.mu <- as.vector(w.beta %*% t(Z))
    Z <- cbind(1, Y.lag, base, aid)
    wo.mu <- as.vector(wo.beta %*% t(Z))
    w <- mean(abs(Y - w.mu))
    wo <- mean(abs(Y - wo.mu))  
    A2 <- rate.taliban.ied.attack.2[t, k] <- (wo-w)*100/w
    
    ####################################################
    Y <- IN$post.B 
    Y.lag <- IN$pre.B 
    XX <- X <- IN$isaf - IN$taliban
    aid <- IN$aid
    base <- IN$base_3km

    #out <- lm(Y ~ X + Y.lag + base + aid + X:base + X:aid)
    #out <- lm(Y ~ X + Y.lag + base + aid + X:Y.lag + X:aid)
    #out <- lm(Y ~ X + Y.lag + base + aid + X:Y.lag + X:base)
    #out <- lm(Y ~ X + Y.lag + base + aid + X:Y.lag)
    out <- lm(Y ~ X + Y.lag + base + aid + X:base)
    #out <- lm(Y ~ X + Y.lag + base + aid + X:aid)
    #out <- lm(Y ~ X + Y.lag + base + aid + X:aid)
    
    w.beta <- coef(out)
    out <- lm(Y ~ Y.lag + base + aid)
    wo.beta <- coef(out)
    Y <- OUT$post.B 
    Y.lag <- OUT$pre.B 
    X<- OUT$isaf - OUT$taliban
    X <- X * (sd(XX)/sd(X))###
    aid <- OUT$aid
    base <- OUT$base_3km
    
    #Z <- cbind(1, X, Y.lag, base, aid, X*base, X*aid)
    #Z <- cbind(1, X, Y.lag, base, aid, X*Y.lag, X*aid)
    #Z <- cbind(1, X, Y.lag, base, aid, X*Y.lag, X*base)
    #Z <- cbind(1, X, Y.lag, base, aid, X*Y.lag)
    Z <- cbind(1, X, Y.lag, base, aid, X*base)
    #Z <- cbind(1, X, Y.lag, base, aid, X*aid)
    #Z <- cbind(1, X, Y.lag, base, aid, X*aid)
    
    w.mu <- w.beta %*% t(Z)
    Z <- cbind(1, Y.lag, base, aid)
    wo.mu <- wo.beta %*% t(Z)
    w <- mean(abs(Y - w.mu))
    wo <- mean(abs(Y - wo.mu))  
    B2 <- rate.taliban.ied.found.2[t, k] <- (wo-w)*100/w

    ####################################################
    Y <- IN$post.C 
    Y.lag <- IN$pre.C 
    XX <- X <- IN$isaf - IN$taliban
    aid <- IN$aid
    base <- IN$base_3km

    #out <- lm(Y ~ X + Y.lag + base + aid + X:base + X:aid)
    #out <- lm(Y ~ X + Y.lag + base + aid + X:Y.lag + X:aid)
    #out <- lm(Y ~ X + Y.lag + base + aid + X:Y.lag + X:base)
    #out <- lm(Y ~ X + Y.lag + base + aid + X:Y.lag)
    out <- lm(Y ~ X + Y.lag + base + aid + X:base)
    #out <- lm(Y ~ X + Y.lag + base + aid + X:aid)
    #out <- lm(Y ~ X + Y.lag + base + aid + X:Y.lag)
    
    w.beta <- coef(out)
    out <- lm(Y ~ Y.lag + base + aid)
    wo.beta <- coef(out)
    Y <- OUT$post.C 
    Y.lag <- OUT$pre.C 
    X<- OUT$isaf - OUT$taliban
    X <- X * (sd(XX)/sd(X))###
    aid <- OUT$aid
    base <- OUT$base_3km
    
    #Z <- cbind(1, X, Y.lag, base, aid, X*base, X*aid)
    #Z <- cbind(1, X, Y.lag, base, aid, X*Y.lag, X*aid)
    #Z <- cbind(1, X, Y.lag, base, aid, X*Y.lag, X*base)
    #Z <- cbind(1, X, Y.lag, base, aid, X*Y.lag)
    Z <- cbind(1, X, Y.lag, base, aid, X*base)
    #Z <- cbind(1, X, Y.lag, base, aid, X*aid)
    #Z <- cbind(1, X, Y.lag, base, aid, X*Y.lag)
    
    w.mu <- w.beta %*% t(Z)
    Z <- cbind(1, Y.lag, base, aid)
    wo.mu <- wo.beta %*% t(Z)
    w <- mean(abs(Y - w.mu))
    wo <- mean(abs(Y - wo.mu))  
    C2 <- rate.taliban.nonied.2[t, k] <- (wo-w)*100/w

########################################################################
    Y <- IN$post.A
    Y.lag <- IN$pre.A
    XX <- X<- IN$isaf - IN$taliban
    aid <- IN$aid
    base <- IN$base_3km
    
    #out <- lm(Y ~ X + Y.lag + base + aid + X:base + X:aid)
    #out <- lm(Y ~ X + Y.lag + base + aid + X:Y.lag + X:aid)
    #out <- lm(Y ~ X + Y.lag + base + aid + X:Y.lag + X:base)
    #out <- lm(Y ~ X + Y.lag + base + aid + X:Y.lag)
    #out <- lm(Y ~ X + Y.lag + base + aid + X:base)
    out <- lm(Y ~ X + Y.lag + base + aid + X:aid)
    #out <- lm(Y ~ X + Y.lag + base + aid)
    
    w.beta <- coef(out)
    out <- lm(Y ~ Y.lag + base + aid)
    wo.beta <- coef(out)
    Y <- OUT$post.A
    Y.lag <- OUT$pre.A
    X<- OUT$isaf - OUT$taliban
    X <- X * (sd(XX)/sd(X))### 
    aid <- OUT$aid
    base <- OUT$base_3km

    #Z <- cbind(1, X, Y.lag, base, aid, X*base, X*aid)
    #Z <- cbind(1, X, Y.lag, base, aid, X*Y.lag, X*aid)
    #Z <- cbind(1, X, Y.lag, base, aid, X*Y.lag, X*base)
    #Z <- cbind(1, X, Y.lag, base, aid, X*Y.lag)
    #Z <- cbind(1, X, Y.lag, base, aid, X*base)
    Z <- cbind(1, X, Y.lag, base, aid, X*aid)
    #Z <- cbind(1, X, Y.lag, base, aid)    
    
    w.mu <- as.vector(w.beta %*% t(Z))
    Z <- cbind(1, Y.lag, base, aid)
    wo.mu <- as.vector(wo.beta %*% t(Z))
    w <- mean(abs(Y - w.mu))
    wo <- mean(abs(Y - wo.mu))  
    A3 <- rate.taliban.ied.attack.3[t, k] <- (wo-w)*100/w
    
    ####################################################
    Y <- IN$post.B 
    Y.lag <- IN$pre.B 
    XX <- X <- IN$isaf - IN$taliban
    aid <- IN$aid
    base <- IN$base_3km

    #out <- lm(Y ~ X + Y.lag + base + aid + X:base + X:aid)
    #out <- lm(Y ~ X + Y.lag + base + aid + X:Y.lag + X:aid)
    #out <- lm(Y ~ X + Y.lag + base + aid + X:Y.lag + X:base)
    #out <- lm(Y ~ X + Y.lag + base + aid + X:Y.lag)
    #out <- lm(Y ~ X + Y.lag + base + aid + X:base)
    out <- lm(Y ~ X + Y.lag + base + aid + X:aid)
    #out <- lm(Y ~ X + Y.lag + base + aid + X:aid)
    
    w.beta <- coef(out)
    out <- lm(Y ~ Y.lag + base + aid)
    wo.beta <- coef(out)
    Y <- OUT$post.B 
    Y.lag <- OUT$pre.B 
    X<- OUT$isaf - OUT$taliban
    X <- X * (sd(XX)/sd(X))###
    aid <- OUT$aid
    base <- OUT$base_3km
    
    #Z <- cbind(1, X, Y.lag, base, aid, X*base, X*aid)
    #Z <- cbind(1, X, Y.lag, base, aid, X*Y.lag, X*aid)
    #Z <- cbind(1, X, Y.lag, base, aid, X*Y.lag, X*base)
    #Z <- cbind(1, X, Y.lag, base, aid, X*Y.lag)
    #Z <- cbind(1, X, Y.lag, base, aid, X*base)
    Z <- cbind(1, X, Y.lag, base, aid, X*aid)
    #Z <- cbind(1, X, Y.lag, base, aid, X*aid)
    
    w.mu <- w.beta %*% t(Z)
    Z <- cbind(1, Y.lag, base, aid)
    wo.mu <- wo.beta %*% t(Z)
    w <- mean(abs(Y - w.mu))
    wo <- mean(abs(Y - wo.mu))  
    B3 <- rate.taliban.ied.found.3[t, k] <- (wo-w)*100/w

    ####################################################
    Y <- IN$post.C 
    Y.lag <- IN$pre.C 
    XX <- X <- IN$isaf - IN$taliban
    aid <- IN$aid
    base <- IN$base_3km

    #out <- lm(Y ~ X + Y.lag + base + aid + X:base + X:aid)
    #out <- lm(Y ~ X + Y.lag + base + aid + X:Y.lag + X:aid)
    #out <- lm(Y ~ X + Y.lag + base + aid + X:Y.lag + X:base)
    #out <- lm(Y ~ X + Y.lag + base + aid + X:Y.lag)
    #out <- lm(Y ~ X + Y.lag + base + aid + X:base)
    out <- lm(Y ~ X + Y.lag + base + aid + X:aid)
    #out <- lm(Y ~ X + Y.lag + base + aid + X:Y.lag)
    
    w.beta <- coef(out)
    out <- lm(Y ~ Y.lag + base + aid)
    wo.beta <- coef(out)
    Y <- OUT$post.C 
    Y.lag <- OUT$pre.C 
    X<- OUT$isaf - OUT$taliban
    X <- X * (sd(XX)/sd(X))###
    aid <- OUT$aid
    base <- OUT$base_3km
    
    #Z <- cbind(1, X, Y.lag, base, aid, X*base, X*aid)
    #Z <- cbind(1, X, Y.lag, base, aid, X*Y.lag, X*aid)
    #Z <- cbind(1, X, Y.lag, base, aid, X*Y.lag, X*base)
    #Z <- cbind(1, X, Y.lag, base, aid, X*Y.lag)
    #Z <- cbind(1, X, Y.lag, base, aid, X*base)
    Z <- cbind(1, X, Y.lag, base, aid, X*aid)
    #Z <- cbind(1, X, Y.lag, base, aid, X*Y.lag)
    
    w.mu <- w.beta %*% t(Z)
    Z <- cbind(1, Y.lag, base, aid)
    wo.mu <- wo.beta %*% t(Z)
    w <- mean(abs(Y - w.mu))
    wo <- mean(abs(Y - wo.mu))  
    C3 <- rate.taliban.nonied.3[t, k] <- (wo-w)*100/w
    
    
    if (verbose!=0&t%%verbose == 0){
      cat("t", t, "k", k, "\n") 
    }
  }
}
##############################################################
png(file = "pred_int_2.png", width = 2200, height = 2200)


ccc <- 3
mmm <- 3
ddd <- 3
fff <- 3
XLAB <- list("Distance window (km)", cex = ccc)
YLAB <- list("Time window (months)", cex = ccc)
COL.MIN <- -36
COL.MAX <- 36
BY <- 5


Q <- rate.taliban.ied.attack.1
MAIN <- list("IED attacks", cex = mmm)
elevation.df <- data.frame(x = rep(1:K, each = T), y = rep(seq(1, MONTH[T], 1), K), z = as.vector(Q))
elevation.loess <- loess(z ~ x * y, data = elevation.df)
elevation.fit <- expand.grid(list(x = seq(1, K, 0.05), y = seq(1, MONTH[T], 0.05)))
z <- predict(elevation.loess, newdata = elevation.fit)
elevation.fit$z <- as.numeric(z)
P1 <- levelplot(z ~ x*y, data = elevation.fit,
                xlab = XLAB, ylab = YLAB,
                main = MAIN,
                col.regions = colorRampPalette(c('dark red','white','dark blue')),
                contour = TRUE,
                at = seq(COL.MIN, COL.MAX, BY), 
                labels = list(TRUE, cex = fff),
                colorkey = FALSE,
                aspect = "xy", scales=list(x=list(cex = ddd),y=list(cex = ddd))
                )


Q <- rate.taliban.ied.found.1
MAIN <- list("IED found", cex = mmm)
elevation.df <- data.frame(x = rep(1:K, each = T), y = rep(seq(1, MONTH[T], 1), K), z = as.vector(Q))
elevation.loess <- loess(z ~ x * y, data = elevation.df)
elevation.fit <- expand.grid(list(x = seq(1, K, 0.05), y = seq(1, MONTH[T], 0.05)))
z <- predict(elevation.loess, newdata = elevation.fit)
elevation.fit$z <- as.numeric(z)
P2 <- levelplot(z ~ x*y, data = elevation.fit,
                xlab = XLAB, ylab = YLAB,
                main = MAIN,
                col.regions = colorRampPalette(c('dark red','white','dark blue')),
                contour = TRUE,
                at = seq(COL.MIN, COL.MAX, BY), 
                labels = list(TRUE, cex = fff),
                colorkey = FALSE,
                aspect = "xy", scales=list(x=list(cex = ddd),y=list(cex = ddd))
                )

Q <- rate.taliban.nonied.1
MAIN <- list("Non-IED attacks", cex = mmm)
elevation.df <- data.frame(x = rep(1:K, each = T), y = rep(seq(1, MONTH[T], 1), K), z = as.vector(Q))
elevation.loess <- loess(z ~ x * y, data = elevation.df)
elevation.fit <- expand.grid(list(x = seq(1, K, 0.05), y = seq(1, MONTH[T], 0.05)))
z <- predict(elevation.loess, newdata = elevation.fit)
elevation.fit$z <- as.numeric(z)
P3 <- levelplot(z ~ x*y, data = elevation.fit,
                xlab = XLAB, ylab = YLAB,
                main = MAIN,
                col.regions = colorRampPalette(c('dark red','white','dark blue')),
                contour = TRUE,
                at = seq(COL.MIN, COL.MAX, BY), 
                labels = list(TRUE, cex = fff),
                colorkey = list(TRUE, labels = list(cex = fff)),
                aspect = "xy", scales=list(x=list(cex = ddd),y=list(cex = ddd))
                )


Q <- rate.taliban.ied.attack.2
MAIN <- list("IED attacks", cex = mmm)
elevation.df <- data.frame(x = rep(1:K, each = T), y = rep(seq(1, MONTH[T], 1), K), z = as.vector(Q))
elevation.loess <- loess(z ~ x * y, data = elevation.df)
elevation.fit <- expand.grid(list(x = seq(1, K, 0.05), y = seq(1, MONTH[T], 0.05)))
z <- predict(elevation.loess, newdata = elevation.fit)
elevation.fit$z <- as.numeric(z)
P4 <- levelplot(z ~ x*y, data = elevation.fit,
                xlab = XLAB, ylab = YLAB,
                main = MAIN,
                col.regions = colorRampPalette(c('dark red','white','dark blue')),
                contour = TRUE,
                at = seq(COL.MIN, COL.MAX, BY), 
                labels = list(TRUE, cex = fff),
                colorkey = FALSE,
                aspect = "xy", scales=list(x=list(cex = ddd),y=list(cex = ddd))
                )


Q <- rate.taliban.ied.found.2
MAIN <- list("IED found", cex = mmm)
elevation.df <- data.frame(x = rep(1:K, each = T), y = rep(seq(1, MONTH[T], 1), K), z = as.vector(Q))
elevation.loess <- loess(z ~ x * y, data = elevation.df)
elevation.fit <- expand.grid(list(x = seq(1, K, 0.05), y = seq(1, MONTH[T], 0.05)))
z <- predict(elevation.loess, newdata = elevation.fit)
elevation.fit$z <- as.numeric(z)
P5 <- levelplot(z ~ x*y, data = elevation.fit,
                xlab = XLAB, ylab = YLAB,
                main = MAIN,
                col.regions = colorRampPalette(c('dark red','white','dark blue')),
                contour = TRUE,
                at = seq(COL.MIN, COL.MAX, BY), 
                labels = list(TRUE, cex = fff),
                colorkey = FALSE,
                aspect = "xy", scales=list(x=list(cex = ddd),y=list(cex = ddd))
                )

Q <- rate.taliban.nonied.2
MAIN <- list("Non-IED attacks", cex = mmm)
elevation.df <- data.frame(x = rep(1:K, each = T), y = rep(seq(1, MONTH[T], 1), K), z = as.vector(Q))
elevation.loess <- loess(z ~ x * y, data = elevation.df)
elevation.fit <- expand.grid(list(x = seq(1, K, 0.05), y = seq(1, MONTH[T], 0.05)))
z <- predict(elevation.loess, newdata = elevation.fit)
elevation.fit$z <- as.numeric(z)
P6 <- levelplot(z ~ x*y, data = elevation.fit,
                xlab = XLAB, ylab = YLAB,
                main = MAIN,
                col.regions = colorRampPalette(c('dark red','white','dark blue')),
                contour = TRUE,
                at = seq(COL.MIN, COL.MAX, BY), 
                labels = list(TRUE, cex = fff),
                colorkey = list(TRUE, labels = list(cex = fff)),
                aspect = "xy", scales=list(x=list(cex = ddd),y=list(cex = ddd))
                )


Q <- rate.taliban.ied.attack.3
MAIN <- list("IED attacks", cex = mmm)
elevation.df <- data.frame(x = rep(1:K, each = T), y = rep(seq(1, MONTH[T], 1), K), z = as.vector(Q))
elevation.loess <- loess(z ~ x * y, data = elevation.df)
elevation.fit <- expand.grid(list(x = seq(1, K, 0.05), y = seq(1, MONTH[T], 0.05)))
z <- predict(elevation.loess, newdata = elevation.fit)
elevation.fit$z <- as.numeric(z)
P7 <- levelplot(z ~ x*y, data = elevation.fit,
                xlab = XLAB, ylab = YLAB,
                main = MAIN,
                col.regions = colorRampPalette(c('dark red','white','dark blue')),
                contour = TRUE,
                at = seq(COL.MIN, COL.MAX, BY), 
                labels = list(TRUE, cex = fff),
                colorkey = FALSE,
                aspect = "xy", scales=list(x=list(cex = ddd),y=list(cex = ddd))
                )


Q <- rate.taliban.ied.found.3
MAIN <- list("IED found", cex = mmm)
elevation.df <- data.frame(x = rep(1:K, each = T), y = rep(seq(1, MONTH[T], 1), K), z = as.vector(Q))
elevation.loess <- loess(z ~ x * y, data = elevation.df)
elevation.fit <- expand.grid(list(x = seq(1, K, 0.05), y = seq(1, MONTH[T], 0.05)))
z <- predict(elevation.loess, newdata = elevation.fit)
elevation.fit$z <- as.numeric(z)
P8 <- levelplot(z ~ x*y, data = elevation.fit,
                xlab = XLAB, ylab = YLAB,
                main = MAIN,
                col.regions = colorRampPalette(c('dark red','white','dark blue')),
                contour = TRUE,
                at = seq(COL.MIN, COL.MAX, BY), 
                labels = list(TRUE, cex = fff),
                colorkey = FALSE,
                aspect = "xy", scales=list(x=list(cex = ddd),y=list(cex = ddd))
                )

Q <- rate.taliban.nonied.3
MAIN <- list("Non-IED attacks", cex = mmm)
elevation.df <- data.frame(x = rep(1:K, each = T), y = rep(seq(1, MONTH[T], 1), K), z = as.vector(Q))
elevation.loess <- loess(z ~ x * y, data = elevation.df)
elevation.fit <- expand.grid(list(x = seq(1, K, 0.05), y = seq(1, MONTH[T], 0.05)))
z <- predict(elevation.loess, newdata = elevation.fit)
elevation.fit$z <- as.numeric(z)
P9 <- levelplot(z ~ x*y, data = elevation.fit,
                xlab = XLAB, ylab = YLAB,
                main = MAIN,
                col.regions = colorRampPalette(c('dark red','white','dark blue')),
                contour = TRUE,
                at = seq(COL.MIN, COL.MAX, BY), 
                labels = list(TRUE, cex = fff),
                colorkey = list(TRUE, labels = list(cex = fff)),
                aspect = "xy", scales=list(x=list(cex = ddd),y=list(cex = ddd))
                )


ad1 <- 0.04
ad2 <- -0.02
ad3 <- 0.0
print(P1, position = c(0, 0.69 + ad2, 0.315, 0.97 -ad1), more = TRUE)
print(P2, position = c(0.315, 0.69  + ad2, 0.63, 0.97 -ad1), more = TRUE)
print(P3, position = c(0.63, 0.69  + ad2, 1 , 0.97 -ad1), more = TRUE)
 
print(P4, position = c(0, 0.35 + ad2, 0.315, 0.64 -ad1), more = TRUE)
print(P5, position = c(0.315, 0.35  + ad2, 0.63 , 0.64 -ad1), more = TRUE)
print(P6, position = c(0.63 , 0.35  + ad2, 1, 0.64 -ad1), more = TRUE)
 
print(P7, position = c(0, 0.02 + ad2, 0.315, 0.30 -ad1), more = TRUE)
print(P8, position = c(0.315, 0.02  + ad2 , 0.63, 0.30 -ad1), more = TRUE)
print(P9, position = c(0.63 , 0.02  + ad2 , 1, 0.30 -ad1))


panel.text(x = 1050, y = 95, labels = expression(paste("Support " %*% " Aid")), cex = 4, font = 3)

panel.text(x = 1050, y = 825, labels = expression(paste("Support " %*% " Bases")), cex = 4, font = 3)

panel.text(x = 1050, y = 1575, labels = expression(paste("Support " %*% " Past attacks")), cex = 4, font = 3)


dev.off()
##############################################################





########################################################################
########################################################################
###  Figure 25: Interaction term II
########################################################################
########################################################################
verbose <- 1
K <- 60
T <- 10
MONTH <- 1:10
########################################################################
rate.taliban.ied.attack.1 <- matrix(NA, T, K)
rate.taliban.ied.found.1 <- matrix(NA, T, K)
rate.taliban.nonied.1 <- matrix(NA, T, K)

rate.taliban.ied.attack.2 <- matrix(NA, T, K)
rate.taliban.ied.found.2 <- matrix(NA, T, K)
rate.taliban.nonied.2 <- matrix(NA, T, K)

rate.taliban.ied.attack.3 <- matrix(NA, T, K)
rate.taliban.ied.found.3 <- matrix(NA, T, K)
rate.taliban.nonied.3 <- matrix(NA, T, K)
########################################################################
for (k in 1:K){
  for (t in 1:T){
    taliban <- village.data$support.taliban
    isaf <- village.data$support.isaf
    pre.A <- pre.ied.attack[[t]][, k]
    pre.B <- pre.ied.found[[t]][, k]
    pre.C <- pre.nonied[[t]][, k]
    post.A <- post.ied.attack[[t]][, k]
    post.B <- post.ied.found[[t]][, k]
    post.C <- post.nonied[[t]][, k]
    data <- data.frame(pre.A, post.A, pre.B, post.B, pre.C, post.C, taliban, isaf, village.data)
    IN <- na.omit(data)
    
    model <- lm(taliban ~
                village.data$log.population +
                village.data$log.altitude +
                village.data$control +
                village.data$pakistan +
                village.data$talibsharia06)
    beta.taliban <- coef(model)

    model <- lm(isaf ~
                village.data$log.population +
                village.data$log.altitude +
                village.data$control +
                village.data$pakistan +
                village.data$talibsharia06)
    beta.isaf <- coef(model)


    pre.A <- other.pre.ied.attack[[t]][, k]
    pre.B <- other.pre.ied.found[[t]][, k]
    pre.C <- other.pre.nonied[[t]][, k]
    post.A <- other.post.ied.attack[[t]][, k]
    post.B <- other.post.ied.found[[t]][, k]
    post.C <- other.post.nonied[[t]][, k]
    
    X <- cbind(1,
           log(out.PPL$population),
           log(out.PPL$altitude),
           out.PPL$control,
           out.PPL$pakistan,
           out.PPL$talibsharia06)
    
    isaf <- as.vector(beta.isaf %*% t(X))
    taliban <- as.vector(beta.taliban %*% t(X))
    
    data <- data.frame(pre.A, post.A, pre.B, post.B, pre.C, post.C, taliban, isaf, out.PPL)
    OUT <- data
    OUT <- OUT[out.PPL$population != 0, ]
########################################################################
    w <- wo <- rep(NA, nrow(OUT))
########################################################################
    Y <- IN$post.A
    Y.lag <- IN$pre.A
    XX <- X<- IN$isaf - IN$taliban
    aid <- IN$aid
    base <- IN$base_3km
    
    out <- lm(Y ~ X + Y.lag + base + aid + X:base + X:aid)
    #out <- lm(Y ~ X + Y.lag + base + aid + X:Y.lag + X:aid)
    #out <- lm(Y ~ X + Y.lag + base + aid + X:Y.lag + X:base)
    #out <- lm(Y ~ X + Y.lag + base + aid + X:Y.lag)
    #out <- lm(Y ~ X + Y.lag + base + aid + X:base)
    #out <- lm(Y ~ X + Y.lag + base + aid + X:aid)
    #out <- lm(Y ~ X + Y.lag + base + aid)
    
    w.beta <- coef(out)
    out <- lm(Y ~ Y.lag + base + aid)
    wo.beta <- coef(out)
    Y <- OUT$post.A
    Y.lag <- OUT$pre.A
    X<- OUT$isaf - OUT$taliban
    X <- X * (sd(XX)/sd(X))### 
    aid <- OUT$aid
    base <- OUT$base_3km

    Z <- cbind(1, X, Y.lag, base, aid, X*base, X*aid)
    #Z <- cbind(1, X, Y.lag, base, aid, X*Y.lag, X*aid)
    #Z <- cbind(1, X, Y.lag, base, aid, X*Y.lag, X*base)
    #Z <- cbind(1, X, Y.lag, base, aid, X*Y.lag)
    #Z <- cbind(1, X, Y.lag, base, aid, X*base)
    #Z <- cbind(1, X, Y.lag, base, aid, X*aid)
    #Z <- cbind(1, X, Y.lag, base, aid)    
    
    w.mu <- as.vector(w.beta %*% t(Z))
    Z <- cbind(1, Y.lag, base, aid)
    wo.mu <- as.vector(wo.beta %*% t(Z))
    w <- mean(abs(Y - w.mu))
    wo <- mean(abs(Y - wo.mu))  
    A1 <- rate.taliban.ied.attack.1[t, k] <- (wo-w)*100/w
    
    ####################################################
    Y <- IN$post.B 
    Y.lag <- IN$pre.B 
    XX <- X <- IN$isaf - IN$taliban
    aid <- IN$aid
    base <- IN$base_3km

    out <- lm(Y ~ X + Y.lag + base + aid + X:base + X:aid)
    #out <- lm(Y ~ X + Y.lag + base + aid + X:Y.lag + X:aid)
    #out <- lm(Y ~ X + Y.lag + base + aid + X:Y.lag + X:base)
    #out <- lm(Y ~ X + Y.lag + base + aid + X:Y.lag)
    #out <- lm(Y ~ X + Y.lag + base + aid + X:base)
    #out <- lm(Y ~ X + Y.lag + base + aid + X:aid)
    #out <- lm(Y ~ X + Y.lag + base + aid + X:aid)
    
    w.beta <- coef(out)
    out <- lm(Y ~ Y.lag + base + aid)
    wo.beta <- coef(out)
    Y <- OUT$post.B 
    Y.lag <- OUT$pre.B 
    X<- OUT$isaf - OUT$taliban
    X <- X * (sd(XX)/sd(X))###
    aid <- OUT$aid
    base <- OUT$base_3km
    
    Z <- cbind(1, X, Y.lag, base, aid, X*base, X*aid)
    #Z <- cbind(1, X, Y.lag, base, aid, X*Y.lag, X*aid)
    #Z <- cbind(1, X, Y.lag, base, aid, X*Y.lag, X*base)
    #Z <- cbind(1, X, Y.lag, base, aid, X*Y.lag)
    #Z <- cbind(1, X, Y.lag, base, aid, X*base)
    #Z <- cbind(1, X, Y.lag, base, aid, X*aid)
    #Z <- cbind(1, X, Y.lag, base, aid, X*aid)
    
    w.mu <- w.beta %*% t(Z)
    Z <- cbind(1, Y.lag, base, aid)
    wo.mu <- wo.beta %*% t(Z)
    w <- mean(abs(Y - w.mu))
    wo <- mean(abs(Y - wo.mu))  
    B1 <- rate.taliban.ied.found.1[t, k] <- (wo-w)*100/w

    ####################################################
    Y <- IN$post.C 
    Y.lag <- IN$pre.C 
    XX <- X <- IN$isaf - IN$taliban
    aid <- IN$aid
    base <- IN$base_3km

    out <- lm(Y ~ X + Y.lag + base + aid + X:base + X:aid)
    #out <- lm(Y ~ X + Y.lag + base + aid + X:Y.lag + X:aid)
    #out <- lm(Y ~ X + Y.lag + base + aid + X:Y.lag + X:base)
    #out <- lm(Y ~ X + Y.lag + base + aid + X:Y.lag)
    #out <- lm(Y ~ X + Y.lag + base + aid + X:base)
    #out <- lm(Y ~ X + Y.lag + base + aid + X:aid)
    #out <- lm(Y ~ X + Y.lag + base + aid + X:Y.lag)
    
    w.beta <- coef(out)
    out <- lm(Y ~ Y.lag + base + aid)
    wo.beta <- coef(out)
    Y <- OUT$post.C 
    Y.lag <- OUT$pre.C 
    X<- OUT$isaf - OUT$taliban
    X <- X * (sd(XX)/sd(X))###
    aid <- OUT$aid
    base <- OUT$base_3km
    
    Z <- cbind(1, X, Y.lag, base, aid, X*base, X*aid)
    #Z <- cbind(1, X, Y.lag, base, aid, X*Y.lag, X*aid)
    #Z <- cbind(1, X, Y.lag, base, aid, X*Y.lag, X*base)
    #Z <- cbind(1, X, Y.lag, base, aid, X*Y.lag)
    #Z <- cbind(1, X, Y.lag, base, aid, X*base)
    #Z <- cbind(1, X, Y.lag, base, aid, X*aid)
    #Z <- cbind(1, X, Y.lag, base, aid, X*Y.lag)
    
    w.mu <- w.beta %*% t(Z)
    Z <- cbind(1, Y.lag, base, aid)
    wo.mu <- wo.beta %*% t(Z)
    w <- mean(abs(Y - w.mu))
    wo <- mean(abs(Y - wo.mu))  
    C1 <- rate.taliban.nonied.1[t, k] <- (wo-w)*100/w


########################################################################
    Y <- IN$post.A
    Y.lag <- IN$pre.A
    XX <- X<- IN$isaf - IN$taliban
    aid <- IN$aid
    base <- IN$base_3km
    
    #out <- lm(Y ~ X + Y.lag + base + aid + X:base + X:aid)
    out <- lm(Y ~ X + Y.lag + base + aid + X:Y.lag + X:aid)
    #out <- lm(Y ~ X + Y.lag + base + aid + X:Y.lag + X:base)
    #out <- lm(Y ~ X + Y.lag + base + aid + X:Y.lag)
    #out <- lm(Y ~ X + Y.lag + base + aid + X:base)
    #out <- lm(Y ~ X + Y.lag + base + aid + X:aid)
    #out <- lm(Y ~ X + Y.lag + base + aid)
    
    w.beta <- coef(out)
    out <- lm(Y ~ Y.lag + base + aid)
    wo.beta <- coef(out)
    Y <- OUT$post.A
    Y.lag <- OUT$pre.A
    X<- OUT$isaf - OUT$taliban
    X <- X * (sd(XX)/sd(X))### 
    aid <- OUT$aid
    base <- OUT$base_3km

    #Z <- cbind(1, X, Y.lag, base, aid, X*base, X*aid)
    Z <- cbind(1, X, Y.lag, base, aid, X*Y.lag, X*aid)
    #Z <- cbind(1, X, Y.lag, base, aid, X*Y.lag, X*base)
    #Z <- cbind(1, X, Y.lag, base, aid, X*Y.lag)
    #Z <- cbind(1, X, Y.lag, base, aid, X*base)
    #Z <- cbind(1, X, Y.lag, base, aid, X*aid)
    #Z <- cbind(1, X, Y.lag, base, aid)    
    
    w.mu <- as.vector(w.beta %*% t(Z))
    Z <- cbind(1, Y.lag, base, aid)
    wo.mu <- as.vector(wo.beta %*% t(Z))
    w <- mean(abs(Y - w.mu))
    wo <- mean(abs(Y - wo.mu))  
    A2 <- rate.taliban.ied.attack.2[t, k] <- (wo-w)*100/w
    
    ####################################################
    Y <- IN$post.B 
    Y.lag <- IN$pre.B 
    XX <- X <- IN$isaf - IN$taliban
    aid <- IN$aid
    base <- IN$base_3km

    #out <- lm(Y ~ X + Y.lag + base + aid + X:base + X:aid)
    out <- lm(Y ~ X + Y.lag + base + aid + X:Y.lag + X:aid)
    #out <- lm(Y ~ X + Y.lag + base + aid + X:Y.lag + X:base)
    #out <- lm(Y ~ X + Y.lag + base + aid + X:Y.lag)
    #out <- lm(Y ~ X + Y.lag + base + aid + X:base)
    #out <- lm(Y ~ X + Y.lag + base + aid + X:aid)
    #out <- lm(Y ~ X + Y.lag + base + aid + X:aid)
    
    w.beta <- coef(out)
    out <- lm(Y ~ Y.lag + base + aid)
    wo.beta <- coef(out)
    Y <- OUT$post.B 
    Y.lag <- OUT$pre.B 
    X<- OUT$isaf - OUT$taliban
    X <- X * (sd(XX)/sd(X))###
    aid <- OUT$aid
    base <- OUT$base_3km
    
    #Z <- cbind(1, X, Y.lag, base, aid, X*base, X*aid)
    Z <- cbind(1, X, Y.lag, base, aid, X*Y.lag, X*aid)
    #Z <- cbind(1, X, Y.lag, base, aid, X*Y.lag, X*base)
    #Z <- cbind(1, X, Y.lag, base, aid, X*Y.lag)
    #Z <- cbind(1, X, Y.lag, base, aid, X*base)
    #Z <- cbind(1, X, Y.lag, base, aid, X*aid)
    #Z <- cbind(1, X, Y.lag, base, aid, X*aid)
    
    w.mu <- w.beta %*% t(Z)
    Z <- cbind(1, Y.lag, base, aid)
    wo.mu <- wo.beta %*% t(Z)
    w <- mean(abs(Y - w.mu))
    wo <- mean(abs(Y - wo.mu))  
    B2 <- rate.taliban.ied.found.2[t, k] <- (wo-w)*100/w

    ####################################################
    Y <- IN$post.C 
    Y.lag <- IN$pre.C 
    XX <- X <- IN$isaf - IN$taliban
    aid <- IN$aid
    base <- IN$base_3km

    #out <- lm(Y ~ X + Y.lag + base + aid + X:base + X:aid)
    out <- lm(Y ~ X + Y.lag + base + aid + X:Y.lag + X:aid)
    #out <- lm(Y ~ X + Y.lag + base + aid + X:Y.lag + X:base)
    #out <- lm(Y ~ X + Y.lag + base + aid + X:Y.lag)
    #out <- lm(Y ~ X + Y.lag + base + aid + X:base)
    #out <- lm(Y ~ X + Y.lag + base + aid + X:aid)
    #out <- lm(Y ~ X + Y.lag + base + aid + X:Y.lag)
    
    w.beta <- coef(out)
    out <- lm(Y ~ Y.lag + base + aid)
    wo.beta <- coef(out)
    Y <- OUT$post.C 
    Y.lag <- OUT$pre.C 
    X<- OUT$isaf - OUT$taliban
    X <- X * (sd(XX)/sd(X))###
    aid <- OUT$aid
    base <- OUT$base_3km
    
    #Z <- cbind(1, X, Y.lag, base, aid, X*base, X*aid)
    Z <- cbind(1, X, Y.lag, base, aid, X*Y.lag, X*aid)
    #Z <- cbind(1, X, Y.lag, base, aid, X*Y.lag, X*base)
    #Z <- cbind(1, X, Y.lag, base, aid, X*Y.lag)
    #Z <- cbind(1, X, Y.lag, base, aid, X*base)
    #Z <- cbind(1, X, Y.lag, base, aid, X*aid)
    #Z <- cbind(1, X, Y.lag, base, aid, X*Y.lag)
    
    w.mu <- w.beta %*% t(Z)
    Z <- cbind(1, Y.lag, base, aid)
    wo.mu <- wo.beta %*% t(Z)
    w <- mean(abs(Y - w.mu))
    wo <- mean(abs(Y - wo.mu))  
    C2 <- rate.taliban.nonied.2[t, k] <- (wo-w)*100/w

########################################################################
    Y <- IN$post.A
    Y.lag <- IN$pre.A
    XX <- X<- IN$isaf - IN$taliban
    aid <- IN$aid
    base <- IN$base_3km
    
    #out <- lm(Y ~ X + Y.lag + base + aid + X:base + X:aid)
    #out <- lm(Y ~ X + Y.lag + base + aid + X:Y.lag + X:aid)
    out <- lm(Y ~ X + Y.lag + base + aid + X:Y.lag + X:base)
    #out <- lm(Y ~ X + Y.lag + base + aid + X:Y.lag)
    #out <- lm(Y ~ X + Y.lag + base + aid + X:base)
    #out <- lm(Y ~ X + Y.lag + base + aid + X:aid)
    #out <- lm(Y ~ X + Y.lag + base + aid)
    
    w.beta <- coef(out)
    out <- lm(Y ~ Y.lag + base + aid)
    wo.beta <- coef(out)
    Y <- OUT$post.A
    Y.lag <- OUT$pre.A
    X<- OUT$isaf - OUT$taliban
    X <- X * (sd(XX)/sd(X))### 
    aid <- OUT$aid
    base <- OUT$base_3km

    #Z <- cbind(1, X, Y.lag, base, aid, X*base, X*aid)
    #Z <- cbind(1, X, Y.lag, base, aid, X*Y.lag, X*aid)
    Z <- cbind(1, X, Y.lag, base, aid, X*Y.lag, X*base)
    #Z <- cbind(1, X, Y.lag, base, aid, X*Y.lag)
    #Z <- cbind(1, X, Y.lag, base, aid, X*base)
    #Z <- cbind(1, X, Y.lag, base, aid, X*aid)
    #Z <- cbind(1, X, Y.lag, base, aid)    
    
    w.mu <- as.vector(w.beta %*% t(Z))
    Z <- cbind(1, Y.lag, base, aid)
    wo.mu <- as.vector(wo.beta %*% t(Z))
    w <- mean(abs(Y - w.mu))
    wo <- mean(abs(Y - wo.mu))  
    A3 <- rate.taliban.ied.attack.3[t, k] <- (wo-w)*100/w
    
    ####################################################
    Y <- IN$post.B 
    Y.lag <- IN$pre.B 
    XX <- X <- IN$isaf - IN$taliban
    aid <- IN$aid
    base <- IN$base_3km

    #out <- lm(Y ~ X + Y.lag + base + aid + X:base + X:aid)
    #out <- lm(Y ~ X + Y.lag + base + aid + X:Y.lag + X:aid)
    out <- lm(Y ~ X + Y.lag + base + aid + X:Y.lag + X:base)
    #out <- lm(Y ~ X + Y.lag + base + aid + X:Y.lag)
    #out <- lm(Y ~ X + Y.lag + base + aid + X:base)
    #out <- lm(Y ~ X + Y.lag + base + aid + X:aid)
    #out <- lm(Y ~ X + Y.lag + base + aid + X:aid)
    
    w.beta <- coef(out)
    out <- lm(Y ~ Y.lag + base + aid)
    wo.beta <- coef(out)
    Y <- OUT$post.B 
    Y.lag <- OUT$pre.B 
    X<- OUT$isaf - OUT$taliban
    X <- X * (sd(XX)/sd(X))###
    aid <- OUT$aid
    base <- OUT$base_3km
    
    #Z <- cbind(1, X, Y.lag, base, aid, X*base, X*aid)
    #Z <- cbind(1, X, Y.lag, base, aid, X*Y.lag, X*aid)
    Z <- cbind(1, X, Y.lag, base, aid, X*Y.lag, X*base)
    #Z <- cbind(1, X, Y.lag, base, aid, X*Y.lag)
    #Z <- cbind(1, X, Y.lag, base, aid, X*base)
    #Z <- cbind(1, X, Y.lag, base, aid, X*aid)
    #Z <- cbind(1, X, Y.lag, base, aid, X*aid)
    
    w.mu <- w.beta %*% t(Z)
    Z <- cbind(1, Y.lag, base, aid)
    wo.mu <- wo.beta %*% t(Z)
    w <- mean(abs(Y - w.mu))
    wo <- mean(abs(Y - wo.mu))  
    B3 <- rate.taliban.ied.found.3[t, k] <- (wo-w)*100/w

    ####################################################
    Y <- IN$post.C 
    Y.lag <- IN$pre.C 
    XX <- X <- IN$isaf - IN$taliban
    aid <- IN$aid
    base <- IN$base_3km

    #out <- lm(Y ~ X + Y.lag + base + aid + X:base + X:aid)
    #out <- lm(Y ~ X + Y.lag + base + aid + X:Y.lag + X:aid)
    out <- lm(Y ~ X + Y.lag + base + aid + X:Y.lag + X:base)
    #out <- lm(Y ~ X + Y.lag + base + aid + X:Y.lag)
    #out <- lm(Y ~ X + Y.lag + base + aid + X:base)
    #out <- lm(Y ~ X + Y.lag + base + aid + X:aid)
    #out <- lm(Y ~ X + Y.lag + base + aid + X:Y.lag)
    
    w.beta <- coef(out)
    out <- lm(Y ~ Y.lag + base + aid)
    wo.beta <- coef(out)
    Y <- OUT$post.C 
    Y.lag <- OUT$pre.C 
    X<- OUT$isaf - OUT$taliban
    X <- X * (sd(XX)/sd(X))###
    aid <- OUT$aid
    base <- OUT$base_3km
    
    #Z <- cbind(1, X, Y.lag, base, aid, X*base, X*aid)
    #Z <- cbind(1, X, Y.lag, base, aid, X*Y.lag, X*aid)
    Z <- cbind(1, X, Y.lag, base, aid, X*Y.lag, X*base)
    #Z <- cbind(1, X, Y.lag, base, aid, X*Y.lag)
    #Z <- cbind(1, X, Y.lag, base, aid, X*base)
    #Z <- cbind(1, X, Y.lag, base, aid, X*aid)
    #Z <- cbind(1, X, Y.lag, base, aid, X*Y.lag)
    
    w.mu <- w.beta %*% t(Z)
    Z <- cbind(1, Y.lag, base, aid)
    wo.mu <- wo.beta %*% t(Z)
    w <- mean(abs(Y - w.mu))
    wo <- mean(abs(Y - wo.mu))  
    C3 <- rate.taliban.nonied.3[t, k] <- (wo-w)*100/w
    
    
    if (verbose!=0&t%%verbose == 0){
      cat("t", t, "k", k, "\n") 
    }
  }
}
##############################################################
png(file = "pred_int_1.png", width = 2200, height = 2200)

ccc <- 3
mmm <- 3
ddd <- 3
fff <- 3
XLAB <- list("Distance window (km)", cex = ccc)
YLAB <- list("Time window (months)", cex = ccc)
COL.MIN <- -36
COL.MAX <- 36
BY <- 5


Q <- rate.taliban.ied.attack.1
MAIN <- list("IED attacks", cex = mmm)
elevation.df <- data.frame(x = rep(1:K, each = T), y = rep(seq(1, MONTH[T], 1), K), z = as.vector(Q))
elevation.loess <- loess(z ~ x * y, data = elevation.df)
elevation.fit <- expand.grid(list(x = seq(1, K, 0.05), y = seq(1, MONTH[T], 0.05)))
z <- predict(elevation.loess, newdata = elevation.fit)
elevation.fit$z <- as.numeric(z)
P1 <- levelplot(z ~ x*y, data = elevation.fit,
                xlab = XLAB, ylab = YLAB,
                main = MAIN,
                col.regions = colorRampPalette(c('dark red','white','dark blue')),
                contour = TRUE,
                at = seq(COL.MIN, COL.MAX, BY), 
                labels = list(TRUE, cex = fff),
                colorkey = FALSE,
                aspect = "xy", scales=list(x=list(cex = ddd),y=list(cex = ddd))
                )


Q <- rate.taliban.ied.found.1
MAIN <- list("IED found", cex = mmm)
elevation.df <- data.frame(x = rep(1:K, each = T), y = rep(seq(1, MONTH[T], 1), K), z = as.vector(Q))
elevation.loess <- loess(z ~ x * y, data = elevation.df)
elevation.fit <- expand.grid(list(x = seq(1, K, 0.05), y = seq(1, MONTH[T], 0.05)))
z <- predict(elevation.loess, newdata = elevation.fit)
elevation.fit$z <- as.numeric(z)
P2 <- levelplot(z ~ x*y, data = elevation.fit,
                xlab = XLAB, ylab = YLAB,
                main = MAIN,
                col.regions = colorRampPalette(c('dark red','white','dark blue')),
                contour = TRUE,
                at = seq(COL.MIN, COL.MAX, BY), 
                labels = list(TRUE, cex = fff),
                colorkey = FALSE,
                aspect = "xy", scales=list(x=list(cex = ddd),y=list(cex = ddd))
                )

Q <- rate.taliban.nonied.1
MAIN <- list("Non-IED attacks", cex = mmm)
elevation.df <- data.frame(x = rep(1:K, each = T), y = rep(seq(1, MONTH[T], 1), K), z = as.vector(Q))
elevation.loess <- loess(z ~ x * y, data = elevation.df)
elevation.fit <- expand.grid(list(x = seq(1, K, 0.05), y = seq(1, MONTH[T], 0.05)))
z <- predict(elevation.loess, newdata = elevation.fit)
elevation.fit$z <- as.numeric(z)
P3 <- levelplot(z ~ x*y, data = elevation.fit,
                xlab = XLAB, ylab = YLAB,
                main = MAIN,
                col.regions = colorRampPalette(c('dark red','white','dark blue')),
                contour = TRUE,
                at = seq(COL.MIN, COL.MAX, BY), 
                labels = list(TRUE, cex = fff),
                colorkey = list(TRUE, labels = list(cex = fff)),
                aspect = "xy", scales=list(x=list(cex = ddd),y=list(cex = ddd))
                )


Q <- rate.taliban.ied.attack.2
MAIN <- list("IED attacks", cex = mmm)
elevation.df <- data.frame(x = rep(1:K, each = T), y = rep(seq(1, MONTH[T], 1), K), z = as.vector(Q))
elevation.loess <- loess(z ~ x * y, data = elevation.df)
elevation.fit <- expand.grid(list(x = seq(1, K, 0.05), y = seq(1, MONTH[T], 0.05)))
z <- predict(elevation.loess, newdata = elevation.fit)
elevation.fit$z <- as.numeric(z)
P4 <- levelplot(z ~ x*y, data = elevation.fit,
                xlab = XLAB, ylab = YLAB,
                main = MAIN,
                col.regions = colorRampPalette(c('dark red','white','dark blue')),
                contour = TRUE,
                at = seq(COL.MIN, COL.MAX, BY), 
                labels = list(TRUE, cex = fff),
                colorkey = FALSE,
                aspect = "xy", scales=list(x=list(cex = ddd),y=list(cex = ddd))
                )


Q <- rate.taliban.ied.found.2
MAIN <- list("IED found", cex = mmm)
elevation.df <- data.frame(x = rep(1:K, each = T), y = rep(seq(1, MONTH[T], 1), K), z = as.vector(Q))
elevation.loess <- loess(z ~ x * y, data = elevation.df)
elevation.fit <- expand.grid(list(x = seq(1, K, 0.05), y = seq(1, MONTH[T], 0.05)))
z <- predict(elevation.loess, newdata = elevation.fit)
elevation.fit$z <- as.numeric(z)
P5 <- levelplot(z ~ x*y, data = elevation.fit,
                xlab = XLAB, ylab = YLAB,
                main = MAIN,
                col.regions = colorRampPalette(c('dark red','white','dark blue')),
                contour = TRUE,
                at = seq(COL.MIN, COL.MAX, BY), 
                labels = list(TRUE, cex = fff),
                colorkey = FALSE,
                aspect = "xy", scales=list(x=list(cex = ddd),y=list(cex = ddd))
                )

Q <- rate.taliban.nonied.2
MAIN <- list("Non-IED attacks", cex = mmm)
elevation.df <- data.frame(x = rep(1:K, each = T), y = rep(seq(1, MONTH[T], 1), K), z = as.vector(Q))
elevation.loess <- loess(z ~ x * y, data = elevation.df)
elevation.fit <- expand.grid(list(x = seq(1, K, 0.05), y = seq(1, MONTH[T], 0.05)))
z <- predict(elevation.loess, newdata = elevation.fit)
elevation.fit$z <- as.numeric(z)
P6 <- levelplot(z ~ x*y, data = elevation.fit,
                xlab = XLAB, ylab = YLAB,
                main = MAIN,
                col.regions = colorRampPalette(c('dark red','white','dark blue')),
                contour = TRUE,
                at = seq(COL.MIN, COL.MAX, BY), 
                labels = list(TRUE, cex = fff),
                colorkey = list(TRUE, labels = list(cex = fff)),
                aspect = "xy", scales=list(x=list(cex = ddd),y=list(cex = ddd))
                )


Q <- rate.taliban.ied.attack.3
MAIN <- list("IED attacks", cex = mmm)
elevation.df <- data.frame(x = rep(1:K, each = T), y = rep(seq(1, MONTH[T], 1), K), z = as.vector(Q))
elevation.loess <- loess(z ~ x * y, data = elevation.df)
elevation.fit <- expand.grid(list(x = seq(1, K, 0.05), y = seq(1, MONTH[T], 0.05)))
z <- predict(elevation.loess, newdata = elevation.fit)
elevation.fit$z <- as.numeric(z)
P7 <- levelplot(z ~ x*y, data = elevation.fit,
                xlab = XLAB, ylab = YLAB,
                main = MAIN,
                col.regions = colorRampPalette(c('dark red','white','dark blue')),
                contour = TRUE,
                at = seq(COL.MIN, COL.MAX, BY), 
                labels = list(TRUE, cex = fff),
                colorkey = FALSE,
                aspect = "xy", scales=list(x=list(cex = ddd),y=list(cex = ddd))
                )


Q <- rate.taliban.ied.found.3
MAIN <- list("IED found", cex = mmm)
elevation.df <- data.frame(x = rep(1:K, each = T), y = rep(seq(1, MONTH[T], 1), K), z = as.vector(Q))
elevation.loess <- loess(z ~ x * y, data = elevation.df)
elevation.fit <- expand.grid(list(x = seq(1, K, 0.05), y = seq(1, MONTH[T], 0.05)))
z <- predict(elevation.loess, newdata = elevation.fit)
elevation.fit$z <- as.numeric(z)
P8 <- levelplot(z ~ x*y, data = elevation.fit,
                xlab = XLAB, ylab = YLAB,
                main = MAIN,
                col.regions = colorRampPalette(c('dark red','white','dark blue')),
                contour = TRUE,
                at = seq(COL.MIN, COL.MAX, BY), 
                labels = list(TRUE, cex = fff),
                colorkey = FALSE,
                aspect = "xy", scales=list(x=list(cex = ddd),y=list(cex = ddd))
                )

Q <- rate.taliban.nonied.3
MAIN <- list("Non-IED attacks", cex = mmm)
elevation.df <- data.frame(x = rep(1:K, each = T), y = rep(seq(1, MONTH[T], 1), K), z = as.vector(Q))
elevation.loess <- loess(z ~ x * y, data = elevation.df)
elevation.fit <- expand.grid(list(x = seq(1, K, 0.05), y = seq(1, MONTH[T], 0.05)))
z <- predict(elevation.loess, newdata = elevation.fit)
elevation.fit$z <- as.numeric(z)
P9 <- levelplot(z ~ x*y, data = elevation.fit,
                xlab = XLAB, ylab = YLAB,
                main = MAIN,
                col.regions = colorRampPalette(c('dark red','white','dark blue')),
                contour = TRUE,
                at = seq(COL.MIN, COL.MAX, BY), 
                labels = list(TRUE, cex = fff),
                colorkey = list(TRUE, labels = list(cex = fff)),
                aspect = "xy", scales=list(x=list(cex = ddd),y=list(cex = ddd))
                )




ad1 <- 0.04
ad2 <- -0.02
ad3 <- 0.0
print(P1, position = c(0, 0.69 + ad2, 0.315, 0.97 -ad1), more = TRUE)
print(P2, position = c(0.315, 0.69  + ad2, 0.63, 0.97 -ad1), more = TRUE)
print(P3, position = c(0.63, 0.69  + ad2, 1 , 0.97 -ad1), more = TRUE)
 
print(P4, position = c(0, 0.35 + ad2, 0.315, 0.64 -ad1), more = TRUE)
print(P5, position = c(0.315, 0.35  + ad2, 0.63 , 0.64 -ad1), more = TRUE)
print(P6, position = c(0.63 , 0.35  + ad2, 1, 0.64 -ad1), more = TRUE)
 
print(P7, position = c(0, 0.02 + ad2, 0.315, 0.30 -ad1), more = TRUE)
print(P8, position = c(0.315, 0.02  + ad2 , 0.63, 0.30 -ad1), more = TRUE)
print(P9, position = c(0.63 , 0.02  + ad2 , 1, 0.30 -ad1))

panel.text(x = 1050, y = 95, labels = expression(paste("Support " %*% " Aid,   Support " %*% " Bases")), cex = 4, font = 3)

panel.text(x = 1050, y = 825, labels = expression(paste("Support " %*% " Past attacks,   Support " %*% " Aid")), cex = 4, font = 3)

panel.text(x = 1050, y = 1575, labels = expression(paste("Support " %*% " Past attacks,   Support " %*% " Bases")), cex = 4, font = 3)


dev.off()
##############################################################




########################################################################
########################################################################
###  Figure 26: Interaction term III
########################################################################
########################################################################
verbose <- 1
K <- 60
T <- 10
MONTH <- 1:10
########################################################################
rate.taliban.ied.attack <- matrix(NA, T, K)
rate.taliban.ied.found <- matrix(NA, T, K)
rate.taliban.nonied <- matrix(NA, T, K)
########################################################################
for (k in 1:K){
  for (t in 1:T){
    taliban <- village.data$support.taliban
    isaf <- village.data$support.isaf
    pre.A <- pre.ied.attack[[t]][, k]
    pre.B <- pre.ied.found[[t]][, k]
    pre.C <- pre.nonied[[t]][, k]
    post.A <- post.ied.attack[[t]][, k]
    post.B <- post.ied.found[[t]][, k]
    post.C <- post.nonied[[t]][, k]
    data <- data.frame(pre.A, post.A, pre.B, post.B, pre.C, post.C, taliban, isaf, village.data)
    IN <- na.omit(data)
    
    model <- lm(taliban ~
                village.data$log.population +
                village.data$log.altitude +
                village.data$control +
                village.data$pakistan +
                village.data$talibsharia06)
    beta.taliban <- coef(model)

    model <- lm(isaf ~
                village.data$log.population +
                village.data$log.altitude +
                village.data$control +
                village.data$pakistan +
                village.data$talibsharia06)
    beta.isaf <- coef(model)


    pre.A <- other.pre.ied.attack[[t]][, k]
    pre.B <- other.pre.ied.found[[t]][, k]
    pre.C <- other.pre.nonied[[t]][, k]
    post.A <- other.post.ied.attack[[t]][, k]
    post.B <- other.post.ied.found[[t]][, k]
    post.C <- other.post.nonied[[t]][, k]
    
    X <- cbind(1,
           log(out.PPL$population),
           log(out.PPL$altitude),
           out.PPL$control,
           out.PPL$pakistan,
           out.PPL$talibsharia06)
    
    isaf <- as.vector(beta.isaf %*% t(X))
    taliban <- as.vector(beta.taliban %*% t(X))
    
    data <- data.frame(pre.A, post.A, pre.B, post.B, pre.C, post.C, taliban, isaf, out.PPL)
    OUT <- data
    OUT <- OUT[out.PPL$population != 0, ]
########################################################################
    w <- wo <- rep(NA, nrow(OUT))
########################################################################
    Y <- IN$post.A
    Y.lag <- IN$pre.A
    XX <- X<- IN$isaf - IN$taliban
    aid <- IN$aid
    base <- IN$base_3km
    out <- lm(Y ~ X + Y.lag + base + aid + X:Y.lag + X:base + X:aid)
    w.beta <- coef(out)
    out <- lm(Y ~ Y.lag + base + aid)
    wo.beta <- coef(out)
    Y <- OUT$post.A
    Y.lag <- OUT$pre.A
    X<- OUT$isaf - OUT$taliban
    X <- X * (sd(XX)/sd(X))### 
    aid <- OUT$aid
    base <- OUT$base_3km
    Z <- cbind(1, X, Y.lag, base, aid, X*Y.lag, X*base, X*aid)
    w.mu <- as.vector(w.beta %*% t(Z))
    Z <- cbind(1, Y.lag, base, aid)
    wo.mu <- as.vector(wo.beta %*% t(Z))
    w <- mean(abs(Y - w.mu))
    wo <- mean(abs(Y - wo.mu))  
    A1 <- rate.taliban.ied.attack[t, k] <- (wo-w)*100/w
    
    Y <- IN$post.B 
    Y.lag <- IN$pre.B 
    XX <- X <- IN$isaf - IN$taliban
    aid <- IN$aid
    base <- IN$base_3km
    out <- lm(Y ~ X + Y.lag + base + aid + X:Y.lag + X:base + X:aid)
    w.beta <- coef(out)
    out <- lm(Y ~ Y.lag + base + aid)
    wo.beta <- coef(out)
    Y <- OUT$post.B 
    Y.lag <- OUT$pre.B 
    X<- OUT$isaf - OUT$taliban
    X <- X * (sd(XX)/sd(X))### 
    aid <- OUT$aid
    base <- OUT$base_3km
    Z <- cbind(1, X, Y.lag, base, aid, X*Y.lag, X*base, X*aid)
    w.mu <- w.beta %*% t(Z)
    Z <- cbind(1, Y.lag, base, aid)
    wo.mu <- wo.beta %*% t(Z)
    w <- mean(abs(Y - w.mu))
    wo <- mean(abs(Y - wo.mu))  
    B1 <- rate.taliban.ied.found[t, k] <- (wo-w)*100/w

    Y <- IN$post.C 
    Y.lag <- IN$pre.C 
    XX <- X <- IN$isaf - IN$taliban
    aid <- IN$aid
    base <- IN$base_3km
    out <- lm(Y ~ X + Y.lag + base + aid + X:Y.lag + X:base + X:aid)
    w.beta <- coef(out)
    out <- lm(Y ~ Y.lag + base + aid)
    wo.beta <- coef(out)
    Y <- OUT$post.C 
    Y.lag <- OUT$pre.C 
    X<- OUT$isaf - OUT$taliban
    X <- X * (sd(XX)/sd(X))### 
    aid <- OUT$aid
    base <- OUT$base_3km
    Z <- cbind(1, X, Y.lag, base, aid, X*Y.lag, X*base, X*aid)
    w.mu <- w.beta %*% t(Z)
    Z <- cbind(1, Y.lag, base, aid)
    wo.mu <- wo.beta %*% t(Z)
    w <- mean(abs(Y - w.mu))
    wo <- mean(abs(Y - wo.mu))  
    C1 <- rate.taliban.nonied[t, k] <- (wo-w)*100/w
    
    
    if (verbose!=0&t%%verbose == 0){
      cat("t", t, "k", k, "A1", A1, "B1", B1, "C1", C1, "\n") 
    }
  }
}
##############################################################
png(file = "pred_int_3.png", width = 2200, height = 800) 


ccc <- 3
mmm <- 3
ddd <- 3
fff <- 3
XLAB <- list("Distance window (km)", cex = ccc)
YLAB <- list("Time window (months)", cex = ccc)
COL.MIN <- -37
COL.MAX <- 37
BY <- 5


Q <- rate.taliban.ied.attack
MAIN <- list("IED attacks", cex = mmm)
elevation.df <- data.frame(x = rep(1:K, each = T), y = rep(seq(1, MONTH[T], 1), K), z = as.vector(Q))
elevation.loess <- loess(z ~ x * y, data = elevation.df)
elevation.fit <- expand.grid(list(x = seq(1, K, 0.05), y = seq(1, MONTH[T], 0.05)))
z <- predict(elevation.loess, newdata = elevation.fit)
elevation.fit$z <- as.numeric(z)
P1 <- levelplot(z ~ x*y, data = elevation.fit,
                xlab = XLAB, ylab = YLAB,
                main = MAIN,
                col.regions = colorRampPalette(c('dark red','white','dark blue')),
                contour = TRUE,
                at = seq(COL.MIN, COL.MAX, BY), 
                labels = list(TRUE, cex = fff),
                colorkey = FALSE,
                aspect = "xy", scales=list(x=list(cex = ddd),y=list(cex = ddd))
                )

Q <- rate.taliban.ied.found
MAIN <- list("IED found", cex = mmm)
elevation.df <- data.frame(x = rep(1:K, each = T), y = rep(seq(1, MONTH[T], 1), K), z = as.vector(Q))
elevation.loess <- loess(z ~ x * y, data = elevation.df)
elevation.fit <- expand.grid(list(x = seq(1, K, 0.05), y = seq(1, MONTH[T], 0.05)))
z <- predict(elevation.loess, newdata = elevation.fit)
elevation.fit$z <- as.numeric(z)
P2 <- levelplot(z ~ x*y, data = elevation.fit,
                xlab = XLAB, ylab = YLAB,
                main = MAIN,
                col.regions = colorRampPalette(c('dark red','white','dark blue')),
                contour = TRUE,
                at = seq(COL.MIN, COL.MAX, BY), 
                labels = list(TRUE, cex = fff),
                colorkey = FALSE,
                aspect = "xy", scales=list(x=list(cex = ddd),y=list(cex = ddd))
                )

Q <- rate.taliban.nonied
MAIN <- list("Non-IED attacks", cex = mmm)
elevation.df <- data.frame(x = rep(1:K, each = T), y = rep(seq(1, MONTH[T], 1), K), z = as.vector(Q))
elevation.loess <- loess(z ~ x * y, data = elevation.df)
elevation.fit <- expand.grid(list(x = seq(1, K, 0.05), y = seq(1, MONTH[T], 0.05)))
z <- predict(elevation.loess, newdata = elevation.fit)
elevation.fit$z <- as.numeric(z)
P3 <- levelplot(z ~ x*y, data = elevation.fit,
                xlab = XLAB, ylab = YLAB,
                main = MAIN,
                col.regions = colorRampPalette(c('dark red','white','dark blue')),
                contour = TRUE,
                at = seq(COL.MIN, COL.MAX, BY), 
                labels = list(TRUE, cex = fff),
                colorkey = list(TRUE, labels = list(cex = fff)),
                aspect = "xy", scales=list(x=list(cex = ddd),y=list(cex = ddd))
                )

print(P1, position = c(0, 0, 0.315, 0.92), more = TRUE)
print(P2, position = c(0.315, 0, 0.63, 0.92), more = TRUE)
print(P3, position = c(0.63, 0, 1, 0.92))

panel.text(x = 1050, y = 20, labels = expression(paste("Support " %*% "Past attacks,  Support " %*% " Aid,  Support " %*% "Bases")), cex = 4, font = 3)

dev.off()
##############################################################









########################################################################
########################################################################
###  Figure 29: Negative Binomial Model
########################################################################
########################################################################
verbose <- 1
K <- 60
T <- 10
MONTH <- 1:10
########################################################################
rate.taliban.ied.attack.civilian <- matrix(NA, T, K)
rate.taliban.nonied.civilian <- matrix(NA, T, K)
rate.taliban.both.civilian <- matrix(NA, T, K)

rate.taliban.ied.attack.civilian.2 <- matrix(NA, T, K)
rate.taliban.nonied.civilian.2 <- matrix(NA, T, K)
rate.taliban.both.civilian.2 <- matrix(NA, T, K)

ewA <- ewB <- ewC <- ewD <- ewE <- ewF <- ewoA <- ewoB <- ewoC <- ewoD <- ewoE <- ewoF <- ewAA <- ewBB <- ewCC <- ewDD <- ewEE <- ewFF <- matrix(NA, T, K)
########################################################################
K.min <- 7
T.min <- 2
for (k in K.min:K){
  for (t in T.min:T){
    taliban <- village.data$support.taliban
    isaf <- village.data$support.isaf

    pre.D <- pre.ied.attack.civilian[[t]][, k]
    pre.E <- pre.nonied.civilian[[t]][, k]
    post.D <- post.ied.attack.civilian[[t]][, k]
    post.E <- post.nonied.civilian[[t]][, k]
    pre.F <- pre.D + pre.E
    post.F <- post.D + post.E

    isaf.taliban <- isaf - taliban
    data <- data.frame(pre.D, post.D, pre.E, post.E, pre.F, post.F, taliban, isaf, isaf.taliban, village.data)
    IN <- na.omit(data)
    
    model <- lm(taliban ~
                village.data$log.population +
                village.data$log.altitude +
                village.data$control +
                village.data$pakistan +
                village.data$talibsharia06)
    beta.taliban <- coef(model)

    model <- lm(isaf ~
                village.data$log.population +
                village.data$log.altitude +
                village.data$control +
                village.data$pakistan +
                village.data$talibsharia06)
    beta.isaf <- coef(model)

    model <- lm(isaf.taliban~
                village.data$log.population +
                village.data$log.altitude +
                village.data$control +
                village.data$pakistan +
                village.data$talibsharia06)
    beta.isaf.taliban <- coef(model)

    pre.D <- other.pre.ied.attack.civilian[[t]][, k]
    pre.E <- other.pre.nonied.civilian[[t]][, k]
    post.D <- other.post.ied.attack.civilian[[t]][, k]
    post.E <- other.post.nonied.civilian[[t]][, k]
    pre.F <- pre.D + pre.E
    post.F <- post.D + post.E
    
    X <- cbind(1,
               log(out.PPL$population),
               log(out.PPL$altitude),
               out.PPL$control,
               out.PPL$pakistan,
               out.PPL$talibsharia06)
    
    isaf <- as.vector(beta.isaf %*% t(X))
    taliban <- as.vector(beta.taliban %*% t(X))
    isaf.taliban <- as.vector(beta.isaf.taliban %*% t(X))
    
    data <- data.frame(pre.D, post.D, pre.E, post.E, pre.F, post.F, taliban, isaf, isaf.taliban, out.PPL)
    OUT <- data
    OUT <- OUT[out.PPL$population != 0, ]
########################################################################
    w <- wo <- rep(NA, nrow(OUT))
    options(warn = 2)
########################################################################
    D1 <- NULL
    Y <- IN$post.D
    Y.lag <- IN$pre.D
    XX <- X<- IN$isaf - IN$taliban
    aid <- IN$aid
    base <- IN$base_3km
    out <- NULL
    try(out <- glm.nb(Y ~ X + Y.lag + base + aid, control = glm.control(maxit = 500)),silent = TRUE)
    if(!is.null(out)){
        w.beta <- coef(out)
        out <- NULL
        try(out <- glm.nb(Y ~ Y.lag + base + aid, control = glm.control(maxit = 500)),silent = TRUE)
        if(!is.null(out)){
            wo.beta <- coef(out)
            YY <- Y <- OUT$post.D
            YY.lag <- Y.lag <- OUT$pre.D
            X<- OUT$isaf - OUT$taliban   
            X <- X * (sd(XX)/sd(X))###
            aid <- OUT$aid
            base <- OUT$base_3km
            Z <- cbind(1, X, Y.lag, base, aid)
            w.mu <- as.vector(exp(w.beta %*% t(Z)))
            Z <- cbind(1, Y.lag, base, aid)
            wo.mu <- as.vector(exp(wo.beta %*% t(Z)))
            w <- mean(abs(Y - w.mu), na.rm = TRUE)
            wo <- mean(abs(Y - wo.mu), na.rm = TRUE)  
            D1 <- rate.taliban.ied.attack.civilian[t, k] <- (wo-w)*100/w
        }
    }

    E1 <- NULL
    Y <- IN$post.E 
    Y.lag <- IN$pre.E 
    XX <- X<- IN$isaf - IN$taliban
    aid <- IN$aid
    base <- IN$base_3km
    out <- NULL
    try(out <- glm.nb(Y ~ X + Y.lag + base + aid, control = glm.control(maxit = 500)),silent = TRUE)
    if(!is.null(out)){
        w.beta <- coef(out)
        out <- NULL
        try(out <- glm.nb(Y ~ Y.lag + base + aid, control = glm.control(maxit = 500)),silent = TRUE)
        if(!is.null(out)){
            wo.beta <- coef(out)
            YY <- Y <- OUT$post.E 
            YY.lag <- Y.lag <- OUT$pre.E 
            X<- OUT$isaf - OUT$taliban   
            X <- X * (sd(XX)/sd(X))###
            aid <- OUT$aid
            base <- OUT$base_3km
            Z <- cbind(1, X, Y.lag, base, aid)
            w.mu <- as.vector(exp(w.beta %*% t(Z)))
            Z <- cbind(1, Y.lag, base, aid)
            wo.mu <- as.vector(exp(wo.beta %*% t(Z)))
            w <- mean(abs(Y - w.mu), na.rm = TRUE)
            wo <- mean(abs(Y - wo.mu), na.rm = TRUE)  
            E1 <- rate.taliban.nonied.civilian[t, k] <- (wo-w)*100/w
        }
    }

    F1 <- NULL
    Y <- IN$post.F 
    Y.lag <- IN$pre.F 
    XX <- X<- IN$isaf - IN$taliban
    aid <- IN$aid
    base <- IN$base_3km
    out <- NULL
    try(out <- glm.nb(Y ~ X + Y.lag + base + aid, control = glm.control(maxit = 500)),silent = TRUE)
    if(!is.null(out)){
        w.beta <- coef(out)
        out <- NULL
        try(out <- glm.nb(Y ~ Y.lag + base + aid, control = glm.control(maxit = 500)),silent = TRUE)
        if(!is.null(out)){
            wo.beta <- coef(out)
            YY <- Y <- OUT$post.F 
            YY.lag <- Y.lag <- OUT$pre.F 
            X<- OUT$isaf - OUT$taliban   
            X <- X * (sd(XX)/sd(X))###
            aid <- OUT$aid
            base <- OUT$base_3km
            Z <- cbind(1, X, Y.lag, base, aid)
            w.mu <- as.vector(exp(w.beta %*% t(Z)))
            Z <- cbind(1, Y.lag, base, aid)
            wo.mu <- as.vector(exp(wo.beta %*% t(Z)))
            w <- mean(abs(Y - w.mu), na.rm = TRUE)
            wo <- mean(abs(Y - wo.mu), na.rm = TRUE)  
            F1 <- rate.taliban.both.civilian[t, k] <- (wo-w)*100/w
        }
    }

    D2 <- NULL
    Y <- IN$post.D
    Y.lag <- IN$pre.D
    aid <- IN$aid
    base <- IN$base_3km
    X<- cbind(IN$log.population, IN$log.altitude, IN$control, IN$pakistan, IN$talibsharia06) 
    out <- NULL
    try(out <- glm.nb(Y ~ X + Y.lag + base + aid, control = glm.control(maxit = 500)),silent = TRUE)
    if(!is.null(out)){
        w.beta <- coef(out)
        out <- NULL
        try(out <- glm.nb(Y ~ Y.lag + base + aid, control = glm.control(maxit = 500)),silent = TRUE)
        if(!is.null(out)){
            wo.beta <- coef(out)
            YY <- Y <- OUT$post.D
            YY.lag <- Y.lag <- OUT$pre.D
            aid <- OUT$aid
            base <- OUT$base_3km
            X<- cbind(log(OUT$population), log(OUT$altitude), OUT$control, OUT$pakistan, OUT$talibsharia06) 
            Z <- cbind(1, X, Y.lag, base, aid)
            w.mu <- as.vector(exp(w.beta %*% t(Z)))
            Z <- cbind(1, Y.lag, base, aid)
            wo.mu <- as.vector(exp(wo.beta %*% t(Z)))
            w <- mean(abs(Y - w.mu), na.rm = TRUE)
            wo <- mean(abs(Y - wo.mu), na.rm = TRUE)  
            D2 <- rate.taliban.ied.attack.civilian.2[t, k] <- (wo-w)*100/w
        }
    }

    E2 <- NULL
    Y <- IN$post.E 
    Y.lag <- IN$pre.E 
    aid <- IN$aid
    base <- IN$base_3km
    X<- cbind(IN$log.population, IN$log.altitude, IN$control, IN$pakistan, IN$talibsharia06) 
    out <- NULL
    try(out <- glm.nb(Y ~ X + Y.lag + base + aid, control = glm.control(maxit = 500)),silent = TRUE)
    if(!is.null(out)){
        w.beta <- coef(out)
        out <- NULL
        try(out <- glm.nb(Y ~ Y.lag + base + aid, control = glm.control(maxit = 500)),silent = TRUE)
        if(!is.null(out)){
            wo.beta <- coef(out)
            YY <- Y <- OUT$post.E 
            YY.lag <- Y.lag <- OUT$pre.E
            aid <- OUT$aid
            base <- OUT$base_3km
            X<- cbind(log(OUT$population), log(OUT$altitude), OUT$control, OUT$pakistan, OUT$talibsharia06)
            Z <- cbind(1, X, Y.lag, base, aid)
            w.mu <- as.vector(exp(w.beta %*% t(Z)))
            Z <- cbind(1, Y.lag, base, aid)
            wo.mu <- as.vector(exp(wo.beta %*% t(Z)))
            w <- mean(abs(Y - w.mu), na.rm = TRUE)
            wo <- mean(abs(Y - wo.mu), na.rm = TRUE)  
            E2 <- rate.taliban.nonied.civilian.2[t, k] <- (wo-w)*100/w
        }
    }

    F2 <- NULL
    Y <- IN$post.F 
    Y.lag <- IN$pre.F 
    aid <- IN$aid
    base <- IN$base_3km
    X<- cbind(IN$log.population, IN$log.altitude, IN$control, IN$pakistan, IN$talibsharia06) 
    out <- NULL
    try(out <- glm.nb(Y ~ X + Y.lag + base + aid, control = glm.control(maxit = 500)),silent = TRUE)
    if(!is.null(out)){
        w.beta <- coef(out)
        out <- NULL
        try(out <- glm.nb(Y ~ Y.lag + base + aid, control = glm.control(maxit = 500)),silent = TRUE)
        if(!is.null(out)){
            wo.beta <- coef(out)
            YY <- Y <- OUT$post.F
            YY.lag <- Y.lag <- OUT$pre.F 
            aid <- OUT$aid
            base <- OUT$base_3km
            X<- cbind(log(OUT$population), log(OUT$altitude), OUT$control, OUT$pakistan, OUT$talibsharia06)
            Z <- cbind(1, X, Y.lag, base, aid)
            w.mu <- as.vector(exp(w.beta %*% t(Z)))
            Z <- cbind(1, Y.lag, base, aid)
            wo.mu <- as.vector(exp(wo.beta %*% t(Z)))
            w <- mean(abs(Y - w.mu), na.rm = TRUE)
            wo <- mean(abs(Y - wo.mu), na.rm = TRUE)  
            F2 <- rate.taliban.both.civilian.2[t, k] <- (wo-w)*100/w
        }
    }
    
    if (verbose!=0&t%%verbose == 0){
      cat("t", t, "k", k, "D1", D1, "E1", E1, "F1", F1, "D2", D2, "E2", E2, "F2", F2, "\n") 
    }
  }
}
##############################################################
##############################################################
png(file = "pred_support_vs_covariate_civilian_nb.png", width = 1800, height = 1500)

ccc <- 3
mmm <- 3
ddd <- 3
fff <- 3
XLAB <- list("Distance window (km)", cex = ccc)
YLAB <- list("Time window (months)", cex = ccc)
COL.MIN <- -190
COL.MAX <- 190


BY <- 2
BY2 <- 2

Q <- rate.taliban.ied.attack.civilian
MAIN <- list("IED attacks", cex = mmm)
Q <- Q[T.min:T, K.min:K]
elevation.df <- data.frame(x = rep(K.min:K, each = (T - T.min + 1)), y = rep(seq(MONTH[T.min], MONTH[T], 1), (K - K.min + 1)), z = as.vector(Q))
elevation.loess <- loess(z ~ x * y, data = elevation.df)
elevation.fit <- expand.grid(list(x = seq(K.min, K, 0.05), y = seq(MONTH[T.min], MONTH[T], 0.05)))
z <- predict(elevation.loess, newdata = elevation.fit)
elevation.fit$z <- as.numeric(z)
P1 <- levelplot(z ~ x*y, data = elevation.fit,
                xlab = XLAB, ylab = YLAB,
                main = MAIN,
                col.regions = colorRampPalette(c('dark red','white','dark blue')),
                contour = TRUE,
                at = seq(COL.MIN, COL.MAX, BY), 
                labels = list(TRUE, cex = fff),
                colorkey = FALSE,
                #colorkey = list(TRUE, labels = list(cex = fff)),
                aspect = "xy", scales=list(x=list(cex = ddd),y=list(cex = ddd))
                )

BY <- 20
BY2 <- 20

Q <- rate.taliban.nonied.civilian
MAIN <- list("Non-IED attacks", cex = mmm)
Q <- Q[T.min:T, K.min:K]
elevation.df <- data.frame(x = rep(K.min:K, each = (T - T.min + 1)), y = rep(seq(MONTH[T.min], MONTH[T], 1), (K - K.min + 1)), z = as.vector(Q))
elevation.loess <- loess(z ~ x * y, data = elevation.df)
elevation.fit <- expand.grid(list(x = seq(K.min, K, 0.05), y = seq(MONTH[T.min], MONTH[T], 0.05)))
z <- predict(elevation.loess, newdata = elevation.fit)
elevation.fit$z <- as.numeric(z)
P2 <- levelplot(z ~ x*y, data = elevation.fit,
                xlab = XLAB, ylab = YLAB,
                main = MAIN,
                col.regions = colorRampPalette(c('dark red','white','dark blue')),
                contour = TRUE,
                at = seq(COL.MIN, COL.MAX, BY), 
                labels = list(TRUE, cex = fff),
                colorkey = FALSE,
                #colorkey = list(TRUE, labels = list(cex = fff)),
                aspect = "xy", scales=list(x=list(cex = ddd),y=list(cex = ddd))
                )

BY <- 5
BY2 <- 5

Q <- rate.taliban.both.civilian
MAIN <- list("Both", cex = mmm)
Q <- Q[T.min:T, K.min:K]
elevation.df <- data.frame(x = rep(K.min:K, each = (T - T.min + 1)), y = rep(seq(MONTH[T.min], MONTH[T], 1), (K - K.min + 1)), z = as.vector(Q))
elevation.loess <- loess(z ~ x * y, data = elevation.df)
elevation.fit <- expand.grid(list(x = seq(K.min, K, 0.05), y = seq(MONTH[T.min], MONTH[T], 0.05)))
z <- predict(elevation.loess, newdata = elevation.fit)
elevation.fit$z <- as.numeric(z)
P10 <- levelplot(z ~ x*y, data = elevation.fit,
                xlab = XLAB, ylab = YLAB,
                main = MAIN,
                col.regions = colorRampPalette(c('dark red','white','dark blue')),
                contour = TRUE,
                at = seq(COL.MIN, COL.MAX, BY), 
                labels = list(TRUE, cex = fff),
                colorkey = list(TRUE, labels = list(cex = fff)),
                aspect = "xy", scales=list(x=list(cex = ddd),y=list(cex = ddd))
                )




Q <- rate.taliban.ied.attack.civilian.2
MAIN <- list("IED attacks", cex = mmm)
Q <- Q[T.min:T, K.min:K]
elevation.df <- data.frame(x = rep(K.min:K, each = (T - T.min + 1)), y = rep(seq(MONTH[T.min], MONTH[T], 1), (K - K.min + 1)), z = as.vector(Q))
elevation.loess <- loess(z ~ x * y, data = elevation.df)
elevation.fit <- expand.grid(list(x = seq(K.min, K, 0.05), y = seq(MONTH[T.min], MONTH[T], 0.05)))
z <- predict(elevation.loess, newdata = elevation.fit)
elevation.fit$z <- as.numeric(z)
P3 <- levelplot(z ~ x*y, data = elevation.fit,
                xlab = XLAB, ylab = YLAB,
                main = MAIN,
                col.regions = colorRampPalette(c('dark red','white','dark blue')),
                contour = TRUE,
                at = seq(COL.MIN, COL.MAX, BY), 
                labels = list(TRUE, cex = fff),
                colorkey = FALSE,
                #colorkey = list(TRUE, labels = list(cex = fff)),
                aspect = "xy", scales=list(x=list(cex = ddd),y=list(cex = ddd))
                )

BY <- 10
BY2 <- 10

Q <- rate.taliban.nonied.civilian.2
MAIN <- list("Non-IED attacks", cex = mmm)
Q <- Q[T.min:T, K.min:K]
elevation.df <- data.frame(x = rep(K.min:K, each = (T - T.min + 1)), y = rep(seq(MONTH[T.min], MONTH[T], 1), (K - K.min + 1)), z = as.vector(Q))
elevation.loess <- loess(z ~ x * y, data = elevation.df)
elevation.fit <- expand.grid(list(x = seq(K.min, K, 0.05), y = seq(MONTH[T.min], MONTH[T], 0.05)))
z <- predict(elevation.loess, newdata = elevation.fit)
elevation.fit$z <- as.numeric(z)
P4 <- levelplot(z ~ x*y, data = elevation.fit,
                xlab = XLAB, ylab = YLAB,
                main = MAIN,
                col.regions = colorRampPalette(c('dark red','white','dark blue')),
                contour = TRUE,
                at = seq(COL.MIN, COL.MAX, BY), 
                labels = list(TRUE, cex = fff),
                colorkey = FALSE,
                #colorkey = list(TRUE, labels = list(cex = fff)),
                aspect = "xy", scales=list(x=list(cex = ddd),y=list(cex = ddd))
                )

BY <- 5
BY2 <- 5

Q <- rate.taliban.both.civilian.2
MAIN <- list("Both", cex = mmm)
Q <- Q[T.min:T, K.min:K]
elevation.df <- data.frame(x = rep(K.min:K, each = (T - T.min + 1)), y = rep(seq(MONTH[T.min], MONTH[T], 1), (K - K.min + 1)), z = as.vector(Q))
elevation.loess <- loess(z ~ x * y, data = elevation.df)
elevation.fit <- expand.grid(list(x = seq(K.min, K, 0.05), y = seq(MONTH[T.min], MONTH[T], 0.05)))
z <- predict(elevation.loess, newdata = elevation.fit)
elevation.fit$z <- as.numeric(z)
P11 <- levelplot(z ~ x*y, data = elevation.fit,
                xlab = XLAB, ylab = YLAB,
                main = MAIN,
                col.regions = colorRampPalette(c('dark red','white','dark blue')),
                contour = TRUE,
                at = seq(COL.MIN, COL.MAX, BY), 
                labels = list(TRUE, cex = fff),
                colorkey = list(TRUE, labels = list(cex = fff)),
                aspect = "xy", scales=list(x=list(cex = ddd),y=list(cex = ddd))
                )


ad1 <- 0.02
print(P1, position = c(0, 0.5 + ad1, 0.315, 0.95), more = TRUE)
print(P2, position = c(0.315, 0.5 + ad1, 0.63, 0.95), more = TRUE)
print(P10, position = c(0.63, 0.5 + ad1, 1, 0.95), more = TRUE)
print(P3, position = c(0, 0 + ad1, 0.315, 0.45), more = TRUE)
print(P4, position = c(0.315, 0 + ad1, 0.63, 0.45), more = TRUE)
print(P11, position = c(0.63, 0 + ad1, 1, 0.45))

panel.text(x = 850, y = 40, labels = "Prediction improvement due to support measure", cex = 4, font = 1)
panel.text(x = 850, y = 800, labels = "Prediction improvement due to covariates", cex = 4, font = 1)


dev.off()
##############################################################










